{ "cells": [ { "cell_type": "markdown", "id": "5e1ba1bf", "metadata": {}, "source": [ "# 1-D Heat Equation" ] }, { "cell_type": "markdown", "id": "b0ab5e8b", "metadata": {}, "source": [ "**강좌**: *기초 전산유체역학*" ] }, { "cell_type": "markdown", "id": "f1c5b57e", "metadata": {}, "source": [ "## Heat Equation (1-D)\n", "열전도에 의한 열유속 (Heat Flux)는 Fourier Law에 의해 표현할 수 있다.\n", "\n", "$$\n", "\\dot{q} = -k \\frac{dT}{dx}\n", "$$\n", "\n", "여기서 $T$는 온도, $k$ 는 열전도 계수이고, $\\dot{q}$ 는 열유속 ($W/m^2$) 이다.\n", "\n", "단위 길이 $\\Delta x$ 에 대해서 내부에너지 변화는 열유속의 구배와 같다.\n", "\n", "$$\n", "\\rho C_p \\frac{dT}{dt} = \\lim_{\\Delta x \\rightarrow 0} \n", "\\frac{\\dot{q} + \\dot{q}_x \\Delta x - \\dot{q}}{\\Delta x} = -k \\frac{d^2T}{dx^2}\n", "$$\n", "\n", "이를 정리하면 1차원 Heat equation 이다.\n", "\n", "$$\n", "\\frac{dT}{dt} = \\alpha \\frac{d^2T}{dx^2}.\n", "$$\n", "\n", ":::{figure-md} markdown-fig\n", "\"conduction-fig\"\n", "\n", "Heat Conduction (image from wikimedia.org)\n", ":::" ] }, { "cell_type": "markdown", "id": "e1393359", "metadata": {}, "source": [ "## Finite Difference Method\n", "Linear Wave 방정식 해석과 같이 Finite Difference 를 이용해서 이 편미분 방정식을 차분한다.\n", "\n", "### FTCS (Forward difference in Time, Central difference in Space)\n", "$$\n", "\\frac{T(x,t+\\Delta t) - T(x, t)}{\\Delta t} \n", "=\n", "\\alpha \\frac {T(x + \\Delta x, t) -2 T(x, t) + T(x - \\Delta x, t)}{\\Delta x^2} + O((\\Delta t), (\\Delta x)^2)\n", "$$\n", "\n", "여기서 시간과 공간을 각각 M, N 개로 나누고, $(n, i)$ 점에서 값을 다음과 같이 간단히 표시하자\n", "\n", "$$\n", "T(x_i, t_n) = T_i^n\n", "$$\n", "\n", ":::{figure-md} markdown-fig\n", "\"stencil-fig\"\n", "\n", "FTCS Stencil\n", ":::\n", "\n", "그 결과 위 차분식 (FDE, Finite Difference Equation)은 다음과 같이 표현할 수 있다.\n", "\n", "$$\n", "\\frac{T_i^{n+1}- T_i^n}{\\Delta t} \n", "=\n", "\\alpha \\frac {T_{i+1}^n -2 T_i^n + T_{i-1}^n}{\\Delta x^2}\n", "$$" ] }, { "cell_type": "markdown", "id": "4e589657", "metadata": {}, "source": [ "#### 계산 격자 구성, Solution array 구성\n", "계산 영역을 $n_x + 1$개의 점으로 나누어보자.\n", "\n", "즉 격자점은 다음과 같다.\n", "\n", "$$\n", "x_j = 0, \\Delta x, -1 + 2 \\Delta x, ..., 1 - \\Delta x, 1 ~~(1 \\le j \\le n_x +1)\n", "$$\n", "\n", ":::{figure-md} Grid\n", "\n", "\n", "Grid\n", ":::\n", "\n", "첫번째 격자점과 마지막 격자점은 경계조건으로 부여할 수 있다.\n", "\n", "$$\n", "T_0^n = T_0 \\\\\n", "T_{n_x+1}^n = T_N\n", "$$\n", "\n", "경계 조건을 위해서 Solution array는 경계조건을 포함해서 격자점 개수와 같이 $n_x+1$ 으로 구성한다. \n", "\n", "계산은 양 끝점을 제외한 $T_1, T_2,... T_{n_x}$ 에 대해서 수행한다.\n", "\n", "$$\n", "T_i^{n+1}\n", "=\n", "T_i^n \n", "+\n", "\\frac{\\alpha \\Delta t}{\\Delta x^2}\n", "\\left (\n", "{T_{i+1}^n -2 T_i^n + T_{i-1}^n}\n", "\\right )\n", "$$\n", "\n" ] }, { "cell_type": "markdown", "id": "94f3a24d", "metadata": {}, "source": [ "### 예제\n", "$\\alpha=0.1$ 이고 $x\\in [0, 1]$ 에 대해서 해석한다.\n", "\n", "- 초기 조건 : $T(x,0) = 100$\n", "- 경계 조건 : $T(0, t) = 100$, $T(1, t) = 300$\n", "\n", "$\\Delta t =0.01, \\Delta x=0.1$ 에 $t=1$ 일때 온도 분포를 구하시오.\n", "\n", "*참고* Exact Solution은 변수 분리법에 따라 다음과 같다.\n", "\n", "$$\n", "T = T_0 + \\Delta T x + \\sum_{n=1}^{\\infty} \\frac{2 \\Delta T (-1)^n}{n \\pi} e^{-n^2 \\pi^2 \\alpha t} \\sin(n \\pi x)\n", "$$" ] }, { "cell_type": "code", "execution_count": 1, "id": "22d56f9d", "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "from matplotlib import pyplot as plt\n", "\n", "import numpy as np\n", "\n", "plt.style.use('ggplot')\n", "plt.rcParams['figure.dpi'] = 150" ] }, { "cell_type": "code", "execution_count": 2, "id": "69f5a01c", "metadata": {}, "outputs": [], "source": [ "def cs(nx, u, du, adtdx2):\n", " \"\"\"\n", " Central difference in space\n", " \n", " Parameters\n", " ----------\n", " nx : integer\n", " number of elements in solution array\n", " u : array\n", " solution\n", " du : array\n", " difference of current and next solution\n", " adtdx2 : float\n", " constant\n", " \"\"\"\n", " # DIY\n", " # du_i = a*dt/dx^2*(u_{i-1} - 2*u_i + u_{i+1}) for i in 1, 2..., nx-2" ] }, { "cell_type": "code", "execution_count": 3, "id": "f6fa300d", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Error: 0.16054\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzwAAAJqCAYAAADuVkWwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAABcSAAAXEgFnn9JSAACWhklEQVR4nOzdd3xV9f3H8fe9yc2eZA9C2Al7D5WhKIggINbRVvqz1bZa6yjDhYuKigNcFasoLhwgoCIgCoKyQVaYhr1XAoHs5I7z+yPNhZjByLg3yev5ePRR7/eM+znhm+S+c875HJNhGIYAAAAAoA4yu7oAAAAAAKguBB4AAAAAdRaBBwAAAECdReABAAAAUGcReAAAAADUWQQeAAAAAHUWgQcAAABAnUXgAQAAAFBnEXgAAAAA1FkEHgAAAAB1FoEHAAAAQJ1F4AEAAABQZxF4AAAAANRZBB4AAAAAdZanqwtwV8ePH5dhGC57//DwcElSenq6y2qAe2JuoCLMD5SHuYHyMDdQEVfPD5PJpOjo6Ertg8BTDsMwXBp4zq8DKAtzAxVhfqA8zA2Uh7mBitTm+cElbQAAAADqLAIPAAAAgDqLwAMAAACgziLwAAAAAKizCDwAAAAA6iwCDwAAAIA6i7bUVayqWvYV76c2twBE9agLc8NkMrm6BAAAUE8QeKqAw+FQXl6e8vPzq+xDaEZGhnPfwPnqwtwwmUzy8fGRr6+vzGZONAMAgOpD4KkkwzB09uxZeXp6Kjg4WB4eHlWyX0/Pon8am81WJftD3VEX5obdbldeXp7Onj2rkJAQzvgAAIBqQ+CppLy8PJnNZgUEBFTph7biffFBEL9VF+aGp6enAgIClJmZqby8PPn5+bm6JAAAUEdxLUklFRYWysfHp1Z/+ARcofiytsLCQleXAgAA6jACTyXZbDZZLBZXlwHUShaLpVZfmgcAANwfgacSihsUcHYHuDzF3zu1ueMcAABwbwQeAAAAAHUWgQcAAABAnVVlXdrmzp2rX3/9VQcPHtTZs2dltVoVEhKiVq1aaejQoWrYsGGZ2/38889asGCBDh8+LE9PT7Vo0ULDhw9Xy5Yty32v1NRUzZ49Wzt37pTNZlN8fLwGDBigvn37VtXhAAAAAPWKY9EcmTr2lCksosL1jFNpMjaukvnaITVUWeVUWeD56quvlJ+fr0aNGikhIUGSdOjQIS1dulQrV67UmDFj1LFjxxLbfPTRR5o3b568vLzUrl07Wa1Wbd68WSkpKRo5cqS6detW6n3Wrl2rSZMmyTAMJScnKzAwUFu3btXkyZN14MAB/d///V9VHRIAAABQLzgWzZEx/T0Zi+fKPOq5ckOPcSpNjoljpbTjcki1IvRUWeAZM2aMmjRpIi8vrxLjP/zwg9577z3997//1dtvv+18qvrWrVs1b948BQYGavz48YqJiZEk7dy5U88884wmT56sVq1aKSAgwLmv7OxsTZ48WQ6HQ6NGjVL37t0lSWfOnNFTTz2lefPmqXPnzmrTpk1VHZZL2X74Wkb77rU2Zefm5mratGn64YcftGvXLp09e1Z+fn5q2rSpevfurT/84Q+Ki4tzdZn1wqFDh9SjRw/17NlTM2fOdHU5AADAzZg69pSxeG5RkJk4tszQc37YUUS0TB17uqjaS1Nl9/AkJSWVCjuS1L9/f0VHRysjI0NHjx51jn/77beSpOHDhzvDjiS1aNFC1113nXJzc7VkyZIS+1q8eLFyc3PVpUsXZ9iRpJCQEN1xxx2Sii6tqwtsP3wt++fvyjFxrIxTaeWuVzzxjOnvybFoTg1WWLH169frqquu0rhx47Rp0ya1bNlSgwYNUufOnXXgwAG99tpr6tWrl5YuXerqUmvcypUrFRcXp4ceeqjK9jl9+nTFxcVp4sSJVbZPAABQf5jCImQe9ZwUEe0MPed/Bv1t2KnoLJC7qZGmBcVndTw9i04oFRYWauvWrZKkHj16lFq/eGz9+vUlxotfl7VNp06dZLFYtGXLljrxIEOPzleUO+GKuWvK3r59u2699VadOHFC9913nzZv3qwZM2borbfe0ieffKJNmzZpypQpiomJ0bFjx1xdbr0QHR2tn3/+Wa+//rqrSwEAAG6qrNBjTzsu+//+uzaGHakGAs/PP/+so0ePKiYmRpGRkZKko0ePymq1KigoSGFhYaW2ady4sSTpwIEDJcYPHjwoSWrSpEmpbTw9PZWQkCCr1VriTFJtZQqLlNcjE2pdyjYMQw888IDy8/M1atQoPf744/Lz8yuxjtls1g033KDvvvtO7du3d1Gl9YvFYlGzZs24hBAAAFTo/NCzL9ekB9/4Vqn3/8XtPnNeiiq7h6fYnDlzdOjQIRUUFOjIkSM6dOiQQkND9eCDDzrP9KSnp0tSmWFHknx8fOTv76+cnBzl5eXJ19dXubm5ysnJkSQ1aNCgzO0aNGigPXv2KD09XYmJiRXWOXLkyFJjXl5emjBhgiQpPDz8gsdqGIYyMjLk6elZ5Q8fNZlMMoVHyeuRCSp88VFn6PF6pKi+wvPCjtcjE2QKi6zS979cixcv1o4dOxQbG6uRI0c6z+qVpUGDBiX+LXNzc/Xf//5XX3/9tQ4cOCCLxaLWrVvrzjvv1E033VRq+y5duujQoUM6ceKEpk6dqg8//FAHDhxQRESE7rzzTt13330ymUzavHmzXnzxRa1bt05Wq1W9evXS+PHjS3UOfOCBBzR9+nTNnj1b+fn5eu2117R161ZZLBZdeeWVevzxx9W8efMS27z88st65ZVX9Prrr+v222+vsMbz30OSvvzyS3355ZfOdUePHq0xY8ZIkhYuXKi5c+dq/fr1OnbsmOx2uxo3bqyhQ4fqH//4h7y9vZ1f25tuukkrV66UJE2aNEmTJk1y7rO4roMHD6pr16664oor9NVXX5Wq88svv9THH3+s7du3y263KzExUTfddJP+/ve/y8fHp9yvk9ls1ssvv6xNmzbJZDKpR48eeuqppyrssng+wzBkNpsVERHBA3yrSPG8iIioXb+MUP2YGygPcwOlRETozKMv65Xpm3XML0KjA/+pUb9+rr7PvyyPiGhXV3fJqjzwpKSkaMuWLc7XYWFhuv/++0uclcnPz5ekMu/5Kebt7a2cnBzl5+fL19fXuU3xsvK2OX//rmQYhpSXU4k9/O/Dn6+fLA88Jevr46S04yp85K6icbtdCo+S5YGnJF8/GbnZla656P38K/XBc+HChZKkG2+8scKw81vZ2dkaPny4UlJSFBYW5ryPa/ny5Vq9erXWr1+v8ePHl7ntk08+qY8//lidOnVSQkKCVq1apWeffVa5ubnq27evbr31ViUkJOiqq67Stm3btGDBAv3666/66aef5OvrW2p/3377rT788EO1b99e/fv3144dOzR//nwtX75cX3/9tVq3bn15XxxJ3bp108mTJ7VkyRIlJiaWuBft/P3+61//Um5urlq2bKnk5GRlZWVp48aNeuGFF7R8+XLNmDHD+QeEq6++WjabTWvXrlXr1q1LNO0oPltakdGjR+uTTz6Rj4+PrrrqKvn6+mrlypV6/vnn9cMPP2jmzJllfp1++OEHTZkyRUlJSbr66qu1Y8cOLVq0SBs2bNDPP//sPKMLAABqF4fDoedn/6JjfvGSpLMWf1nl4eKqLl+VB54nn3xSkpSTk6ODBw9q5syZeuaZZ3T77bdr+PDhkv4XBqQKP1gXr1Ndzv8reFnS09MvWINhGHI4HLLZbKWOxcjNluPBP1S6zlLs9vOKPCHrk/+o0t2bX/9MJr+AC69YjuKw27p1a9lstovebvz48UpJSVGvXr30/vvvy9/fX5K0e/du3XzzzZoyZYp69eqlfv36Obcp/veZM2eO5s+f7zyrsHv3bvXv31+TJ0/WjBkz9Mgjj+ivf/2rpKL7x+644w6tWLFCs2fP1m233ebcn8PhkCR98MEHeumll/THP/7R+T4vvPCC3nrrLT344INasGBBqW3sdnuZx1tcY/Gy22+/XQkJCVqyZIm6du1aah4WrzdhwgT17t27xOWA2dnZuu+++7Ro0SLNnDnT+f30j3/8Q2FhYVq7dq0GDBigUaNGldqn/X/zxjCMEnXOmzdPn3zyiaKjozVz5kxnQMrKytKf/vQnrV27Vi+++KKeeOKJUsf87rvv6s0339SwYcOcX4N77rlH8+fP1/vvv+88W1WR4u+htLQ0zvBUkeK/0Kalld/sBPUTcwPlYW7gt76duUjLLfHO17ccXKwOp1OV9vi9NX5Jm8lkKtHg7HJU2z08/v7+Sk5O1mOPPaYmTZpo+vTp2r17tyQ5/1pcUFBQ7vbFjQeKL6c5/7Ka8rYrHv/tJTioORkZGZLKv1yxLLm5ufr8889lNpv1/PPPO8OOJDVr1kwPPvigJGnq1Kllbv/www+XuISqWbNm6tevn/Ly8hQXF+cMO1LRWcW7775bkrRq1aoy99elSxdn2JGKvtHGjBmj2NhYbdmyRevWrbvoY7tc119/fal7nwICAvTMM89IUonQVRnFX9MxY8aUOBsUGBio559/XiaTSZ988kmZjUCGDRvmDDuS5OHhoQceeECStGbNmiqpDwAA1KzUVev0Yd65y9Y6FBzTvWP+esFmWu6s2psWeHp66oorrpBhGM4ua8X3x5w6darMbfLz85WTkyN/f39nOPLz83N+ADx9+nSZ2xWPX8z9N6gel3NmbvPmzcrPz1eHDh3KbEhx8803S5J++eWXMvffq1evUmPFD78ta1mjRo0kSSdPniyznqFDh5Yas1gsuuGGG5x11IS9e/fqvffe0xNPPKGRI0fqoYce0muvveZcVllWq1UbNmyQyWQqEVyKJScnKzk5WdnZ2dq+fXup5X369Ck1VvzvV3zPEgAAqD0y9+3Vy9utspmLLgILLczSv//UV96x8RW2rHZ3VX5JW1kCAwMlSZmZmZKk2NhYWSwWZWZm6tSpU6XOBuzbt0/SuQ+txRo1aqQdO3Zo7969io+PL7HMZrPp4MGDslgsio2Nra5DuXi+/jK//tllb158/4v15HEZ/xkvnTophUXK9M+iS4t+O2ZqUEUhz9f/wutUoLhxRHlhtizFH45/20SgWHBwsIKCgpSZmamsrCwFBQWVWF7Wac7icFzRsvLal/92bhUrru/48eNlLq8qhmHo3//+t6ZMmVJugMzOrvw9WxkZGSosLFRkZGS5Z0UbNmyo7du3lxlgyvraFp+dqwut4QEAqE/sacf1+tzNSmuQJEkyG3aNbuen8Jiie3KLu7cVdwku7+Gk7qhGnsNT/NfhqKgoSUWXFRXfWL169epS6xePde7cucR4p06dyt1mw4YNslqtatOmTYXNEGqKyWSSyS+gUv9TXq6Myc8XBZuIaJnHvCBzfGLR/8a8UJSyT50sWicvr9LvZ/ILqPR9FMU33p/fuOJSvmaXs05N3ftxOWeviu93uRRz5szRu+++q+joaL377rtav3699u/fryNHjjj/GFCV3P3rDgAAqt/XX/+sdf8LO5L0R/90telcslHTb5/TY2ws+/YAd1MlgWfHjh1auXKl88boYjabTd99952WLl0qLy8vXXHFFc5lgwYNkiTNnj27xMMnd+7cqUWLFsnX11fXXHNNif3169dPvr6+WrduXYl7BM6ePatp06ZJkgYPHlwVh+RyxqmTznbUZfU8v9DTcF2luKnA3LlzL7ppQXEQLn7O0m9lZmYqMzNTfn5+Cgi4/IYKF+vw4cNljh85ckRS0UM8i1ksFklF9yH9lt1uv6wbQL/77jtJ0gsvvKBBgwYpOjra+T7lfY0uR2hoqLy8vHTy5Enl5eWVuU7x14KOawAA1F3b1m/VNN9zXV47Fx7VTUOuKnPd4s+gptvulvnaITVVYqVUySVtJ06c0OTJkxUYGKgmTZooMDBQWVlZOnjwoDIyMmSxWPSPf/yjxL017dq10w033KD58+fr4YcfVtu2bWW327V582Y5HA7df//9pT7cBgQE6N5779Wrr76qSZMmqVWrVgoMDNSWLVuUk5OjgQMHqm3btlVxSC5nX7/ygg94+u2pRWPjKplcPPGuvvpqtWzZUqmpqXrjjTfKfN5RsaysLB09elTt2rWTj4+PNm3apL1795a6j2f27NmSilo618RZhTlz5ugvf/lLiTGbzab58+dLKmpqUKw4rJV1T82KFStktVpLjReHl9/+gaDY2bNnJanMSzO//fbbMrcpPqt5KZ3xLBaLOnXqpNWrV+ubb74p9RyhX3/9Vdu3b1dAQIBatWp10fsFAAC1x5m0U3plc54cXkW3oIQXZurBYZ3k4VF+G2pTWITLP3Neiio5w9OqVSvddNNNio2N1cGDB7Vq1Sr9+uuvCggI0PXXX6+JEyeWOLtT7M4779Q//vEPxcXFacuWLUpNTVWbNm00btw49ejRo8z36tGjh8aNG6f27dtr//792rhxo6KionTvvffqz3/+c1Ucjlvw7D9MHr//2wWvjXS3lG0ymfTGG2/Ix8dHEydO1AsvvFDq7IdhGPrhhx80cOBAbdq0SX5+frr99tvlcDg0duzYEuvv2bNHr7/+uiTV2L/vL7/8oi+++KJEva+88oqOHDmiVq1aqWvXrs5lxfN09uzZOnTokHP8wIEDJVo5n6/4DNGePXvKXF4c+D799NMSl9GtWbNGb7/9dpnbFAev8vZZnuKv6SuvvKIDBw44x7Ozs/XEE0/IMAzdcccdbnGZKAAAqFo2m02vfpui0/8LO54Omx7uHKzg0BDXFlbFquQMT2RkpH7/+99f1rZ9+/ZV3759L2mbpKQkPf7445f1frWJZ/9hF/UXe3dL2W3atNEXX3yhv/71r/rPf/6j999/X507d1ZERIQyMzO1efNmpaWlycfHx3kW47HHHtOGDRu0dOlS9ezZUz169FBubq5Wrlyp/Px83XXXXbr22mtrpP4//elPGj16tKZNm+ZslJGamqqAgAC9+uqrJdZt1KiRfve732nmzJnq37+/unfvrtzcXG3YsEH9+vVTQUFBqUvkGjZsqOTkZKWkpGjQoEFq0aKFPDw81L9/f/Xv319/+ctfNGPGDH300UdatWqVkpOTdfz4ca1du1Z///vf9d///rdUzZ06dVJ4eLjmzZun3/3ud0pISJDZbNZtt91WIqD91uDBg/XHP/5Rn376qa655hpdeeWV8vX11apVq3Tq1Cl16tRJo0ePrpovLAAAcCszv1qqTd7nrij5U0iGWrYr3eG2tquRpgWof7p27aoVK1boySefVIcOHbRjxw59++23WrduneLj4zVy5EgtW7bM2TY6ICBAs2bN0ujRo9WgQQMtXLhQa9euVbt27fTWW2/p3//+d43VfuONN+qDDz6Q2WzW999/r2PHjmnAgAH69ttvnc02zvfyyy/rn//8pwICAvTzzz/ryJEjuv/++zV58uRy32PKlCm6/vrrdeDAAc2cOVOff/65s9FD06ZNNX/+fF133XU6ffq0fvjhB+Xk5OjFF190Ptj3t3x8fPTxxx+rd+/e2rZtm2bMmKHPP//8otpXv/TSS3r99dfVpk0brV69WosWLVJ4eLgeeeQRzZgxw9kaHgAA1B2bVqfoi4Jz9yX3tB7VjYOudGFF1cdkXE7rqXrg2LFjF+zKZRiGs612Vd9bUtyW+lLuyUDlPPTQQ/ryyy/15ZdflnkJpruoS3OjOr+H6iuemI7yMDdQHuZG/XPq6En964eDOmspul8+uvCMJt7cRgFBpZtDuXp+mEymMh+FcSk4wwMAAADUE7ZCqyYu2OYMOxaHVQ/3iCgz7NQVBB4AAACgnvhs1s/a5n3ujMnd4ZlqmtzUhRVVPwIPAAAAUA/8smydZjnina9724+q/4CeLqyoZhB4gP957bXXdOTIEbe+fwcAAOBynDx4VK/tPffRP67gtO69uYfM5rofB+r+EQIAAAD1WGF+oV5euEvZnn6SJC97oR7pFSc/fz8XV1YzCDwAAABAHfbxrJ+10yfK+fqemDw1at7IhRXVLAIPAAAAUEetXLxG36qh8/W1Oqp+13V3YUU1j8ADAAAA1EFH9x7Um4e8na8bFZzSX4fXv3uVCTwAAABAHVOQm6eXlhxQrqePJMnXXqCHr0mUj6+PiyureQQeAAAAoI55f+Yy7fOJcL6+L8Gm+MQ4F1bkOgQeAAAAoA756fsV+t4jwfn6BvMx9erb2YUVuRaBBwAAAKgjDu3co7dPBDhfNy1M05+HX+XCilyPwAMAAADUAflZ2Xpx2THlexQ1KvC35enh65rLy9vi4spcy9PVBaBuiYu78LWht9xyi1577bXqLwYAAKCeMAxD78xaqUM+5y5le7CZSdHx0S6syj0QeFAtbrnllnKXdevWrQYruXQTJ07UpEmTNGnSJN12222uLgcAAOCCFs5fpsWWc2FnmOW4ul/Z13UFuRECD6oFZ3AAAABqxr6tOzXlVKjkUfQ6qfCk7ritft+3cz7u4QEAAABqqZwzmXpxbboKPYru0wmy5mr0wFayeHJeoxiBBy41fvx4xcXF6Z577im1LD09XR06dFBCQoLWrVvnHD9x4oQmT56sm2++WZ07d1ZiYqI6dOigu+++W5s2bSr3vXJzc/Xmm29qwIABatGihZo3b66+ffvqqaee0uHDhyVJ3bt316RJkyRJI0eOVFxcnPN/K1eurNqDBwAAqASHw6G3vlqjY94NJEkmw6F/tfJWRHS4iytzL0Q/uNTDDz+sZcuW6dtvv1W/fv1K3PszatQopaWlaeTIkerSpYtz/Pvvv9dzzz2nxMREJSUlKSAgQPv379d3332nRYsW6aOPPlKfPn1KvM+JEyd0++23a+fOnQoJCdGVV14pT09P7d+/X1OnTlXr1q112223adCgQVq2bJm2b9+url27KjEx0bmPyMjIav96AAAAXKzv5izVCq+Gzte/801Tp259KtiifiLwVBPDMJRjdVz29p4OkyTJZrNXVUkXxd9ilslkqrH38/Ly0ltvvaUBAwboySefVI8ePdSwYUN99NFHWrRokTp16qQHH3ywxDZdu3bVwoUL1apVqxLjP/30k/785z/r8ccf1/Lly0scxwMPPKCdO3dq2LBhevnll+Xn5+dctnfvXjkcRf9WTz31lCZOnKjt27fr97//PU0LAACAW9q1cbumZoU7r9dqYz2h27lvp0wEnmqSY3Xoj1/ucnUZl+zTW5orwMuj0vupqD31+++/r+uvv975ulmzZnryySc1duxYPfDAA5owYYKeffZZ+fv7680335Tnb65BTU5OLnO/ffv21eDBgzV79mz9+uuvzvU2btyo5cuXKzIyslTYkaQmTZpc7mECAADUuKxTGXppY6Zs3iGSpBBrtkYNbidPz8p/hquLCDyoFhW1pS4rDN15551avHixfvzxRw0dOlR5eXmaOHFiiUvKzldQUKCffvpJGzdu1OnTp1VYWChJ2rFjhyRp3759zsCzbNkySdJNN91UKuwAAADUJg67TW98s04nvYsuZTMbDo1uH6AG4aEursx9EXhQLS6nLfXEiRPVo0cPZWVl6brrrtPtt99e5no7duzQn//8Zx06dKjcfWVnZzv/++jRo5KkRo0aXXJNAAAA7uSbr37WWu9z9+38IfCU2nbs5cKK3B+Bp5r4W8z69Jbml7198WVcNputqkq6KP4W1zXu++GHH5Sfny9J2r17t3Jzc0udkTEMQ/fcc48OHTqkESNGaMSIEWrUqJH8/f1lMpn0wgsv6D//+Y8Mwyi1/5q8NwkAAKCqbV+boo/zo6X/faTpZD2u4Tf2dm1RtQCBp5qYTKZK3QtTfA2mzVz6g3tdtHfvXj3zzDPy8/NT7969tWDBAj3zzDN66aWXSqy3e/du7d69W+3bt9eECRNK7efgwYOlxmJjYyVJ+/fvr5baAQAAqtuZE2l6eVuBHF7ekqTwwkw9NKyTPMw8ZeZC+ArB5Ww2m+6//37l5uZq3LhxevPNN9W4cWN9+umnWrBgQYl1z5w5I0mKiYkptZ8zZ85o6dKlpcZ79So6zfvVV18pLy/vgvVYLEUP7rLba7ZDHgAAQFnsNptenbtZp72CJEkeDrtGd22g4NAgF1dWOxB44HITJ07Upk2bNGDAAP3hD3+Qn5+fszvbmDFjdPLkSee6jRs3ltls1ooVK7R3717neH5+vh599FFnIDpfx44ddcUVV+jkyZN6+OGHS4Weffv2affu3c7XUVFRkqQ9e/ZU8ZECAABcupmzlmiTz7mmT39qcFbJbZq5sKLahUvaUC0eeuihcpfFxcVpzJgxkqS1a9fqrbfeUkREhF5++WXnOh07dtRDDz2kV155RSNHjtQnn3wik8mk8PBw/f73v9enn36q6667TldeeaV8fHy0du1a2e123XrrrZoxY0ap93zjjTd06623avbs2VqyZIm6devmfPDo9u3bNXHiRDVrVvSDo0+fPvLx8dGUKVOUmpqqqKgomUwm3XPPPc51AAAAasLmFRv0hTXOed9Od/txDRnIfTuXgsCDavHll1+Wu6xVq1YaM2aMsrKy9MADD8hut2vixIkKCwsrsd4DDzygJUuWaMmSJfrggw/0l7/8RZL0wgsvqGnTpvriiy+0YsUKBQYGqlevXnrkkUc0ffr0Mt8zJiZG8+fP15QpUzRv3jz9/PPP8vT0VGxsrO6++25dddW5B3VFR0dr6tSpevXVV7V27Vrl5ORIkoYPH07gAQAANebUkWOauNMhh1fRRVlRhWd1/++6ysx9O5fEZJTVzgo6duxYmZ2+zmcYhk6dOqWwsLAq7wDmqi5tcH91aW5U5/dQfRURESFJSktLc3ElcDfMDZSHueGebAWFemracm3zKWq+5Omw6cVuAWrWMrFG63D1/DCZTGXeu30piIcAAACAm/l85k/OsCNJd0Xm1HjYqSsIPAAAAIAbWffTGs1UgvN1L+O4Bvbv5sKKajcCDwAAAOAmTu4/pNf2n7vNPq4wQ/8Y3oNLvyuBwAMAAAC4AWt+vl5ZtEdZFn9JkpfdqjG9G8rPz8fFldVuBB4AAADADXz85c9K9Y12vv57fIEaN413YUV1A4EHAAAAcLGVC1dojrmR8/U1phO69pouLqyo7iDwAAAAAC50bPc+vXnUz/m6UeFp/X14DxdWVLcQeCqh+OYxHmUEXJ7i7x1uxAQA1FcFObl66eeDyvX0lST52Av08DVN5OPj7eLK6g4CTyV5enrKarW6ugygVrJarc4HqQIAUN8YhqGpXy7VXp8o59g/GxuKbxRdwVa4VASeSvLy8lJ+fj5neYBLZBiG8vPz5eXl5epSAABwiaXfLdMCS6Lz9fWeJ9SrVweX1VNX8afVSvL19VVBQYGys7Pl6+srDw+PKtlvcYAiSOG36sLcsNvtysvLk8PhkK+vr6vLAQCgxh3esUuT04Kdn8abFqbrrluucG1RdRSBp5JMJpOCg4OVl5ens2fPVtmHULO56OSbw+Gokv2h7qgLc8NkMsnHx0fBwcHcvwMAqHfyMzP14srjyveJkCT52fI1ZkBLeXlZXFxZ3UTgqQJms1n+/v7y9/evssATEVH0DZCWllYl+0PdURfmBiEHAFBfGYahd2au1EGfROfYAy09FBMb4bqi6jgCTxWrqg9yxfvhgyF+i7kBAEDttejbn7TYO9H5eoh3mnr26OWyeuqDSgeegoICpaSkaP369dqzZ4/S0tLkcDgUHR2t7t27a/DgwfLx8Smxza233nrB/bZu3VpPP/208/W2bds0bty4ctdv3ry5nnvuucs/EAAAAKAa7du8Q++eCZP+d8t3S2u6/nQb9+1Ut0oHnuXLl+udd96RJDVs2FDt27dXXl6edu7cqRkzZmjFihV65plnFBwc7NymT58+5e5vw4YNysrKUnJycpnLo6KilJSUVOY4AAAA4I5yM87o5V9Oq9AnTJIUaMvV6Btay+JZNQ2vUL5KBx5PT0/1799fgwYNUkxMjHM8IyNDEyZM0L59+/Thhx/qwQcfdC677777ytxXTk6OVq5cKUnq1avsU3tJSUnlbg8AAAC4G4fdrrdmr9ERn0bOsX+18VNkZKgLq6o/Kv0cnj59+ujuu+8uEXYkKTQ0VHfddZckae3atbLZbBfc16pVq2S1WtW8efNS+wMAAABqowVfL9Hy88LO7/xOqXPn0lcsoXpU64NHGzUq+oe1Wq3Kysq64PrLli2TJPXu3bs6ywIAAABqxK71m/V+brTzdRtrmn4/pKcLK6p/qrVL24kTJyRJHh4eCggIqHDd9PR0/frrr/Lw8NAVV5R/89bx48f12WefKSsrS4GBgUpKSlKHDh2czyYBAAAA3EFWWrpeTsmRzbvo0rUQa45GDmknTw8+t9akag088+fPlyR16NBBFkvFD1JatmyZDMNQx44dFRgYWO56qampSk1NLTGWkJCgUaNGXdJlcCNHjiw15uXlpQkTJkiSwsPDL3pf1cHTs+ifpviZK0Ax5gYqwvxAeZgbKA9zo3o4bFa9MHWhTng3lCSZDYee7h2npJbNXFzZpakL86PaAs+GDRu0ZMkSeXh46Lbbbrvg+he6nM3Pz09DhgxR9+7dncFm//79+vzzz7Vr1y6NHz9eL7/8svz8/KruIAAAAIDL8MnUr7T6f2FHkv4Ukafu3du6sKL6y2QYhlHVOz18+LCefPJJ5eTk6M4779QNN9xQ4fp79+7Vo48+Kn9/f7377rsXPBt0PofDoXHjxmnHjh26/fbbNXz48MqWL0k6duyYquFLc9GKU3RaWprLaoB7Ym6gIswPlIe5gfIwN6rejtUbNHaXt+zmopbTHe0n9eQdV8mjFt6C4er5YTKZKt3MrMq/6qdOndLzzz+vnJwcDR48+IJhRzp3dqdHjx6XFHYkyWw2a+jQoZKklJSUSy8YAAAAqCJnj53QyzuszrDTwJqtfw3tVCvDTl1RpV/5zMxMjR8/Xunp6erbt69GjBhxwW0cDscFn71zIdHRRZ0vzpw5c1nbAwAAAJVltxbqtXmbdcorWJLk4bDr4W5hCg6uuHkXqleVBZ68vDy98MILOnLkiLp166Z77rlHJpPpgttt2bJFGRkZioiIUHJy8mW9d05OjiTJx8fnsrYHAAAAKmvWl4u1wffcfTsjwrOV3KqxCyuCVEWBx2q16qWXXtKePXvUvn17PfTQQxfdJrr4crZevXpdVEAqy5o1ayRJjRszoQAAAFDzNi9dq88dCc7X3YyTGnZ9NxdWhGKVDjwOh0Ovv/66tm3bpuTkZI0ePdrZvu5CCgoKtHbtWkkXvpxt4cKFpR5eahiGFi5cqHnz5slkMql///6XdxAAAADAZTp96LAm7jHJYSr6aB1pzdQDw7pc9h/zUbUq3ZZ6wYIFztASGBio9957r8z1RowYoaCgoBJjv/zyi/Lz89W0aVPFxcVV+D5ff/21pk6dqvj4eGe3iIMHD+rkyZMymUy688471aRJk8oeDgAAAHDRbAX5mvh9qs74Fn2W9XTY9PCVMQoM4FEp7qLSgSc7O9v538XBpyy33HJLqcBz/uVsFzJ48GClpKTo8OHD2rJli+x2u0JDQ9WrVy8NHDhQzZrVroc4AQAAoPb74ssl2up77raKv0Tnq3nzNi6sCL9VLc/hqQt4Dg/cFXMDFWF+oDzMDZSHuXH51i9eoWePhsr436VsV+mkRv/h8u9Ld0eunh9u+RweAAAAoK47uXe/Xj3o4ww7sdYzum949zoVduoKAg8AAABwCay5uXpl8V5lWfwlSV4Oqx7umyA/X28XV4ayEHgAAACAi2QYhj6ZsUSpvrHOsb81tKlxYmwFW8GVCDwAAADARVr9w1J9Y2nqfH21OU3X9unguoJwQQQeAAAA4CIc37lbbx4/13W4oTVDfx/eg/t23ByBBwAAALiAwqxsvbT0sHI8fSVJPvZCPdyviXy9LS6uDBdC4AEAAAAqYBiGps78WXt8o51j9zY1KaFhlAurwsUi8AAAAAAVWDZvib7zOnffzgBLuvpe2daFFeFSEHgAAACAchze+qveOhXmfN3Yelp33dTDhRXhUhF4AAAAUC85Fs2RcSqt3OX5Z8/qpdUnle9Z9HwdP3uBHr6+pbwtnjVVIqoAgQcAAAD1jmPRHBnT35Nj4tgyQ4/hcOjdmSt0wDfSOfZPy27FRoeVWhfujcADAACAesfUsacUES2lHS8z9Pz4zY/60aeJ8/Xg46t1Rf/eNV0mqgCBBwAAAPWOKSxC5lHPlRl69m/aoneyznVga551SH8aMVCmsAhXlYtKIPAAAACgXior9OTs2aWX1p1VoYeXJCnAmqMx1zaTdyQtqGsrAg8AAADqrfNDj5F2XG9/s05HfMOdyx9qZlJUs8YurBCVReABAABAvVYcer6P7allUR2d48O9T6jrVZ1cWBmqAoEHAAAA9d72dVs0tdmNztetzuzVH3onubAiVBUCDwAAAOq19E2b9NKRANnMRc/XCSnI1Mjtn8n82pMVPqcHtQOBBwAAAPVWwcF9emnlCZ3xDpIkeTjseridnxoE+5Xbshq1C4EHAAAA9ZIj/bje+3KZUoMbOcfujrOqdZc25basRu1D4AEAAEC9Y5xK0w///UQ/RHdzjvWznNLAq9tLqvg5PahdCDwAAACod36dO19TGg10vm5uO62/39RdJpPJOfbb0GNsXOWKUlFJnq4uAAAAAKhJp/cd0EuOZNk8iz4KB9ty9cjgNvK2lP5oXBx6jI2rZL52SE2XiirAGR4AAADUG4XZ2Xpp0W6d9jqvSUG3MEWEBZW7jSksgrBTixF4AAAAUC8YDoemTv9JO/zinGN/iS1Um+RGFWyF2o7AAwAAgHph0ewf9J1PM+frqz1PadA1HVxXEGoEgQcAAAB13s5Vv+idvHNndpraTuue3zQpQN1E4AEAAECdlnHgoCbscMhqtkiSgmy5enRQa/l40b+rPiDwAAAAoM6yZmfr5R926pR3sCTJbNg1pmsDRYYHu7gy1BQCDwAAAOokw+HQB9OXaJtfvHPszuhCtWuV6LqiUOMIPAAAAKiTlnz1g+b5NHe+7u15SkP6dXBdQXAJAg8AAADqnN2rftHbOeeaFDS2Zeg+mhTUSwQeAAAA1ClnDx7SCzvsKvQoalIQYMvTozck06SgniLwAAAAoM6w5eTo5QW/Kt07RJJkNhwa0zlE0REhLq0LrkPgAQAAQJ1gOBz66IsftcW/oXPsT1H56tCmsQurgqsReAAAAFAnLP36B83xaeF8fZXHaQ27tqMLK4I7IPAAAACg1tu7er3+kx3rfN3Idkb/HN6VJgUg8AAAAKB2O3vosCZsK1Shh5ckKcCer0cHJsnXy+LiyuAOCDwAAACotWx5OZr43Xad8AmVJJkMh0Z1DFJsZIhrC4PbIPAAAACgVjIcDk37bJFS/BOcY3dE5qtT2yYurAruhsADAACAWmn5Nwv1lU9L5+ue5tO6+TqaFKAkAg8AAABqnX1rN+jNrBjn64a2M3qAJgUoA4EHAAAAtUrWkSOasCVfBf9rUuBnz9dj17eUnzdNClAagQcAAAC1hi0vV5PmbdFxnwaSipoUjGwfqLioUBdXBndF4AEAAECtYBiGPv98oTb4JzrHfh+Rp67tm7quKLg9Ag8AAABqhZVzFmqm97kmBd3Np3VL/04urAi1gWdld1BQUKCUlBStX79ee/bsUVpamhwOh6Kjo9W9e3cNHjxYPj4+JbaZMWOGZs6cWe4+hw4dqj/+8Y9lLktNTdXs2bO1c+dO2Ww2xcfHa8CAAerbt29lDwUAAABu6uC6DXrjTJTz02uc7awevL2LzDQpwAVUOvAsX75c77zzjiSpYcOGat++vfLy8rRz507NmDFDK1as0DPPPKPg4OBS27Zs2VLR0dGlxps0Kbt3+tq1azVp0iQZhqHk5GQFBgZq69atmjx5sg4cOKD/+7//q+zhAAAAwM1kHzmq51PylO8TJknysxfo8QHN5e/t5eLKUBtUOvB4enqqf//+GjRokGJizrUGzMjI0IQJE7Rv3z59+OGHevDBB0tt269fv4s+M5Odna3JkyfL4XBo1KhR6t69uyTpzJkzeuqppzRv3jx17txZbdq0qewhAQAAwE3Y83L16twUHQto7Bx7qJ2/4qMbuLAq1CaVvoenT58+uvvuu0uEHUkKDQ3VXXfdJanozIzNZqvU+yxevFi5ubnq0qWLM+xIUkhIiO644w5J0ty5cyv1HgAAAHAfhmHoiy8Wat15Yef2sBx179DMhVWhtqnWpgWNGjWSJFmtVmVlZVVqX+vXr5ck9ejRo9SyTp06yWKxaMuWLSosLKzU+wAAAMA9rP52oWZ4nWtS0NV8WrcNoEkBLk2lL2mryIkTJyRJHh4eCggIKLV869at2r9/vwoLCxUWFqaOHTuWe//OwYMHJZV9f4+np6cSEhK0Z88eHT16VImJiVV3EAAAAKhxh9Zv0usZkc5Pq7G2TD10K00KcOmqNfDMnz9fktShQwdZLKWffLt06dISr6dPn67u3bvrvvvuK9HZLTc3Vzk5OZKkBg3Kvl6zQYMG2rNnj9LT0y8q8IwcObLUmJeXlyZMmCBJCg8Pv+A+qpOnZ9E/TUREhEvrgPthbqAizA+Uh7mB8rjj3Dh76JBe2JitPN+iz2M+9kK9dGsXNU6IcnFl9Y87zo9LVW2BZ8OGDVqyZIk8PDx02223lVgWHR2tESNGqGPHjgoPD1dOTo527NihadOmac2aNXI4HBozZoxz/fz8fOd/e3t7l/l+xePnrwsAAIDaxZaXq2c++UlH/BKdY2N7RKgJYQeXqVoCz+HDh/Xmm2/KMAyNGDGi1BmX3r17l3jt4+Ojq666Sq1bt9bo0aP1yy+/KDU1VS1btlR1mTRpUoXL09PTZRhGtb3/hRSn6LS0NJfVAPfE3EBFmB8oD3MD5XGnuWEYhqZ/+I3W+CU5x24JzVa7FkluUV995Or5YTKZSjVHu1RV3rTg1KlTev7555WTk6PBgwfrhhtuuOhtQ0NDnW2qU1JSnOPnX95WUFBQ5rbF4799yCkAAABqh7VzF+kLSwvn686mDP3++s4urAh1QZUGnszMTI0fP17p6enq27evRowYccn7KE5wGRkZzjE/Pz/5+flJkk6fPl3mdsXjrr73BgAAAJfu8IYUvXY6Qoap6ONpjC1T/7qpszzMNClA5VRZ4MnLy9MLL7ygI0eOqFu3brrnnntkuowuGtnZ2ZJKn6kpbnG9d+/eUtvYbDYdPHhQFotFsbGxl1E9AAAAXCX32DFN2JCpXM+iz38+9kI9em1TBfp6ubgy1AVVEnisVqteeukl7dmzR+3bt9dDDz0ks/nSd20Yhn755RdJpdtPd+pU1HN99erVpbbbsGGDrFar2rRpIy8vvjEAAABqC0d+vl6fs1GHfM91AXugta8S48JcWBXqkkoHHofDoddff13btm1TcnKyRo8e7WxfV5bMzEz9/PPPslqtJcbz8/M1ZcoU7dq1SyEhIerWrVuJ5f369ZOvr6/WrVunNWvWOMfPnj2radOmSZIGDx5c2cMBAABADTEMQ7M+X6DVAef+0D08JFtXdm7uwqpQ11S6S9uCBQu0du1aSVJgYKDee++9MtcbMWKEgoKClJ+fr7feektTp05VfHy8wsLClJubq3379ikrK0v+/v4aOXJkqfbTAQEBuvfee/Xqq69q0qRJatWqlQIDA7Vlyxbl5ORo4MCBatu2bWUPBwAAADVk/bxF+vS8JgUdTBm6Y2APF1aEuqjSgaf4nhtJzuBTlltuuUVBQUEKDAzU0KFDtWvXLh0/flz79++X2WxWZGSk+vTpo8GDB5f7cNEePXpo3Lhxmj17tnbt2iWbzaa4uDgNGDBAV199dWUPBQAAADXk6KYUTUoPl2EpuuAoypalUbfQpABVz2S48mEzbuzYsWM8hwduibmBijA/UB7mBsrjirmRe+KEHpmTqoN+kZIkb0ehXrw6Vo3j6bbrblz9s8Mtn8MDAAAAlMeRn683v/7FGXYk6Z/JvoQdVBsCDwAAAGqEYRj66ovvtDKgmXNsaHC2enehSQGqD4EHAAAANWLj/EWa5tnS+bqdMvR/N3R2YUWoDwg8AAAAqHbHUzZrYloDOUxFHz8jbNkaPawTTQpQ7Qg8AAAAqFb5J0/qhbWnlW3xlyR5Oax67JpGCvb3vsCWQOUReAAAAFBtHAX5+s/sNdrvF+0c+0dLLzVtGOHCqlCfEHgAAABQLQzD0JwvvtOywHNNCW4MytbV3VpWsBVQtQg8AAAAqBYp3y3WRx4tnK/b6IzuHESTAtQsAg8AAACq3InNWzXxRLAcJg9JUrgtW6OHdpAnTQpQwwg8AAAAqFL5aSc1YXWaMr0CJEkWh02PXp2g0AAfF1eG+ojAAwAAgCrjKMjX27NWa69/jHPsnhYWNU+IdGFVqM8IPAAAAKgShmFo3ufz9VPguft2bgjM1rXdaVIA1yHwAAAAoEps+X6xpnqeCzetdEZ3DaZJAVyLwAMAAIBKO7llq14+dq5JQQNbjh4e2p4mBXA5Ag8AAAAqpSAtTS+uOulsUuDpsOnRPg0VGuDr4soAAg8AAAAqwVFYoP/OXqnd/rHOsXuae6plIk0K4B4IPAAAALgshmFowRfztDjg3H071wdk6boeSS6sCiiJwAMAAIDLsu2HxXrPdC7sJBlndPfgLi6sCCiNwAMAAIBLlr5tu14+GiS7uahJQag9Vw8PbS+LB00K4F4IPAAAALgkhenpmrDimM54BUqSPB12PdIrTmGBNCmA+yHwAAAA4KIZVqvenblcu/zjnGN/beqh5MZRLqwKKB+BBwAAABfFMAx9/8VcLQw815TgOv8sXX8FTQrgvgg8AAAAuCg7Fv6kKWrufN3COKu/39jZhRUBF0bgAQAAwAWd2r5DLx32l83sKUkKsefqkRvbyeLBx0m4N2YoAAAAKlR4Ol0vLjuiDO8gSZKHYdfDV8UqPJgmBXB/BB4AAACUy7Ba9d6Xy5QaEO8cu6uxWa2bRLuwKuDiEXgAAABQroVfzNX3AcnO19f4ZekGmhSgFiHwAAAAoEy/Llyid85rUtDMOKt7buwkk4mHi6L2IPAAAACglIwdv+rFQ37OJgVB9jw9cmM7eXt6uLgy4NIQeAAAAFCC9fQpvbT0kE57B0uSzIZDD18RrUiaFKAWIvAAAADAybBZNXXGT9oe0NA59pdGUttmMS6sCrh8BB4AAIB6wrFojoxTaRWu8+P0uZof2Nr5uq9vlgZflVzBFoB7I/AAAADUA45Fc2RMf0+OiWPLDT27flyi/zqaOV83LkjXvTQpQC1H4AEAAKgHTB17ShHRUtrxMkPPmdRfNWG/r6xmiyQp0JqjR69pJB8LTQpQuxF4AAAA6gFTWITMo54rM/QUnErXy0sOKN0nRFJRk4LR7fwVndiwgj0CtQOBBwAAoJ4oK/TYjh7Wa5O/1NbARs71/hSZqw6deLgo6gYCDwAAQD3y29Aza9zz+iawjXN5L88MDbuuswsrBKoWgQcAAKCeKQ49u4IT9HaLm53jifaz+ufwbjQpQJ1C4AEAAKiHjm5K0XOt71Shh5ckKcCao0d7x9KkAHUOgQcAAKCeObN2tZ7d661MrwBJkofDrjHbpiny3Wcv+JweoLYh8AAAANQjedu36Ln1Z3XML8I5NqaVt9pacsptWQ3UZgQeAACAesK6f5cmLdypnUHnOrLd1dSiIdf3LLdlNVDbEXgAAADqAfvRQ5r65TKtDW/tHOsflKe/DO4uqeLn9AC1GYEHAACgjjPsdn0z43vNj73COdbZkqV7BnUo0ZHtt6HH2LjKFeUCVcrT1QUAAACg+hiGoWWffaWPIq5yjjVTpkYP6ywPc+n208Whx9i4SuZrh9RkqUC14AwPAABAHbZ1zny9riTn6yh7jp4Y1l5+XuW3nzaFRRB2UGdU+gxPQUGBUlJStH79eu3Zs0dpaWlyOByKjo5W9+7dNXjwYPn4+DjXdzgcSk1N1fr167V9+3adPHlSubm5CgsLU9u2bTVs2DBFRkaWep9t27Zp3Lhx5dbRvHlzPffcc5U9HAAAgDrjwM/L9EJGjGyWoo98AfZ8PTWohUL9vV1cGVBzKh14li9frnfeeUeS1LBhQ7Vv3155eXnauXOnZsyYoRUrVuiZZ55RcHCwJOnkyZN6+umnJUkNGjRQixYtZDabtXv3bi1atEgrVqzQY489pqSkpDLfLyoqqsxlUVFRlT0UAACAOuNUyib9e7encnz8JEkWh01je8cpPizQxZUBNavSgcfT01P9+/fXoEGDFBMT4xzPyMjQhAkTtG/fPn344Yd68MEHncvat2+vm266Sa1atXKOWa1WTZkyRT/99JPeeOMNvfHGG/L0LF1eUlKS7rvvvsqWDQAAUGfl7t+rZ9ecVrp/rCTJZBj6V4cgtUqMuMCWQN1T6Xt4+vTpo7vvvrtE2JGk0NBQ3XXXXZKktWvXymazSZKio6M1duzYEmFHkiwWi+6++275+fkpPT1dO3furGxpAAAA9Y41/aRe+u5X7ftf2JGkPzc268q2CS6sCnCdam1a0KhR0UOtrFarsrKyLri+l5eXMzidPn26OksDAACocxw5WXp7xnJtDGriHBsSVqihV7Z0YVWAa1VrW+oTJ05Ikjw8PBQQEHDB9R0Oh9LT0yVJISEhZa5z/PhxffbZZ8rKylJgYKCSkpLUoUMHmc00nAMAAPWXYbVq+ifz9GNwB+dYT98c/XlAJ9cVBbiBag088+fPlyR16NBBFovlguuvWLFCZ8+eVVBQkFq2LPsvEampqUpNTS0xlpCQoFGjRpW6rK4iI0eOLDXm5eWlCRMmSJLCw8Mvel/Vofj+pYgIrrVFScwNVIT5gfIwN+o2w+HQzNen6IvADs6x1p45ev6u/vL2rPiPwswNVKQuzI9qCzwbNmzQkiVL5OHhodtuu+2C66enp+vDDz+UJN16662lApKfn5+GDBmi7t27O4PN/v379fnnn2vXrl0aP368Xn75Zfn5+VX5sQAAALiznz+YpjeMltL/niMarxy9/JdrLhh2gPrAZBiGUdU7PXz4sJ588knl5OTozjvv1A033FDh+vn5+XrmmWe0d+9ede3aVWPGjLno93I4HBo3bpx27Nih22+/XcOHD69s+ZKkY8eOqRq+NBetOEWnpaW5rAa4J+YGKsL8QHmYG3XXnh8W6fFjEcr3LHq2Tog9Ty8OaanokIv7IzBzAxVx9fwwmUyXdBVXWao89p86dUrPP/+8cnJyNHjw4AuGHZvNpokTJ2rv3r1KSkoq0b76YpjNZg0dOlSSlJKSctl1AwAA1DYn167Vs4eDnGHH22HVE/0SLzrsAPVBlQaezMxMjR8/Xunp6erbt69GjBhR4foOh0NvvvmmUlJS1KhRIz3yyCPy8vK65PeNjo6WJJ05c+ZyygYAAKh1snb+qn+n5CnDO0iSZDYcGtM1TM3jQl1cGeBeqizw5OXl6YUXXtCRI0fUrVs33XPPPTKZTBVu895772nVqlWKiYnRE088IX9//8t675ycHEmSj4/PZW0PAABQmxQeP6IXFh/QIb8o59jfW3qra1JsBVsB9VOVBB6r1aqXXnpJe/bsUfv27fXQQw9dsE30Z599pkWLFik8PFxPPvmkgoODL/v916xZI0lq3LjxZe8DAACgNrBnntEbs37RtsBGzrFbou26vmtTF1YFuK9KBx6Hw6HXX39d27ZtU3JyskaPHu1sX1eeuXPn6uuvv1ZISIiefPLJi2oBvXDhwlIPLzUMQwsXLtS8efNkMpnUv3//Sh0LAACAOzMKCvTJtO+1LCTJOdY3ME9/vKaVC6sC3Ful21IvWLBAa9eulSQFBgbqvffeK3O9ESNGKCgoSPv379cnn3wiSYqMjNTs2bPLXL9fv35KSjr3zfz1119r6tSpio+Pd3aLOHjwoE6ePCmTyaQ777xTTZo0KXNfAAAAtZ3hsGvex7P1VXBH51g7z2z9c1DnC95GANRnlQ482dnZzv8uDj5lueWWWxQUFKScnBxnu+edO3dq586dZa7funXrEoFn8ODBSklJ0eHDh7VlyxbZ7XaFhoaqV69eGjhwoJo1a1bZQwEAAHBLhmFozeez9L5ve+dYI2XrkWEdZPEg7AAVqZbn8NQFPIcH7oq5gYowP1Ae5kbt9uu38/VkRrwKPYq62TZw5Oqlm1opIqDyDZuYG6iIq+eHWz6HBwAAAFXnyLLlei490hl2/ByFempA8yoJO0B9QOABAABwU2e2bNazqVKmV4AkydOw69Ero9U4MtDFlQG1B4EHAADADeUf3K/nVpzQMd9z3Wz/2dZf7ZtEurAqoPYh8AAAALgZ2+l0TZy7RTsDGzrH7kiQrm6f6LqigFqKwAMAAOBGHDnZeu/zxVob3Nw5NiC0QL+7qqULqwJqLwIPAACAmzBsVn31yRx9F9LOOdbFJ0d/v74dz9oBLhOBBwAAwA0YhqGln8zSx4GdnGPNTNkaM7SjPMyEHeByEXgAAADcwJZZ3+gNzzbO11GOXD0xtJ18PPm4BlQG30EAAAAudmDRj5qQnSCb2VOSFOAo0FM3tFSov5eLKwNqPwIPAACAC51at07/PuCnHIufJMnLYdMTfeMVH+bv4sqAuoHAAwAA4CK5u3fq2fVZSvcJlSSZDEP/6hyi5IZhLq4MqDsIPAAAAC5gPX5UL/6wS/sCYpxjf2nmqStaxbuwKqDuIfAAAADUMEfmWb09c6U2BTd1jg2JtGlIj+YVbAXgchB4AAAAapBRWKAvps3Tj6HnOrJd4Z+nP1/b2oVVAXUXgQcAAKCGGA67Fn40U9ODzz1rJ9kjR/+6sb3MPFgUqBYEHgAAgBpgGIY2fD5L//Xp4ByLM3L0+LD28vLgIxlQXfjuAgAAqAF75n+nl2wtZTd7SJJCHPl66sbWCvLxdHFlQN1G4AEAAKhmJ1Ys1/gT4cr39JYk+TiseuK6xooO9nFxZUDdR+ABAACoRlnbturZbXZleAdJksyGQ6N7RKh5dLCLKwPqBwIPAABANSk8dEAvLD2sQ/5RzrF7kn3VtXm0C6sC6hcCDwAAQDWwZ6Tr9TkbtC0o0Tl2S5w0oHNj1xUF1EMEHgAAgCpm5OXq408XanlIsnPs6uAC/bFPSxdWBdRPBB4AAIAqZNhsmvfRLH0d2tE51t4rV/cNbCcTz9oBahyBBwAAoIoYhqHV02bo/YDOzrFEU44eHdZeFg/CDuAKBB4AAIAqkvrVN5pkbiOHqegjVpgjT08OaSM/i4eLKwPqLwIPAABAFTiyeLGeOxuvQg8vSZKfo1BPDWyu8ABvF1cG1G8EHgAAgEo6s3Gd/r3XS5leAZIkT8Oux3rFKjE8wMWVASDwAAAAVEL+3l0avyZDx33DnWP3tw9Wu8TwCrYCUFMIPAAAAJfJdvKYJn63XbsCGzrHRiR6qG/beBdWBeB8BB4AAIDL4Mg6q/dmLNXakHPP1rk+zKabr2jmwqoA/BaBBwAA4BIZhQX66uM5+i60vXOsi2+e/ta/Nc/aAdwMgQcAAOASGA6Hln40Qx+HdHWONfPI0Zgh7eVhJuwA7obAAwAAcAm2zJipN3w6Ol9HGbl6Ymg7+XjysQpwR3xnAgAAXKQD332nCfnNZDN7SpICHQV6alCyQn0tLq4MQHkIPAAAABfh1OqV+vfREOVY/CRJXg6bxl6ToPhQXxdXBqAiBB4AAIALyN2xTc+m5CvdJ1SSZDIMjewapuS4UBdXBuBCCDwAAAAVsB45qBeX7Ne+gFjn2F0tvdUzKcaFVQG4WAQeAACAcjgyTmny7DXaFNzUOTYk2tCNXZu4sCoAl4LAAwAAUAYjP1eff/qdFjdo6xy7MrBAf74myYVVAbhUBB4AAIDfMGw2LfzwS80I7eIca2XJ1UOD2srMg0WBWoXAAwAAcB7DMLTh0+l6O+Bc2IlTrh4f2l5eHnx0AmobvmsBAADOs+frb/SS2shh8pAkhTry9dSNrRXo7eHiygBcDgIPAADA/xz/eYmezYhVvqe3JMnHYdUTA5oqOsjbxZUBuFwEHgAAAElZmzbo2VSzzngHSZLMhkNjroxWs8hAF1cGoDIIPAAAoN4r3L9bz686qcP+Uc6xe9sEqEuTCBdWBaAqEHgAAECd5lg0R8aptHKX29NP6LW5m7U9KNE5dmtDs/p3SKiB6gBUN8/K7qCgoEApKSlav3699uzZo7S0NDkcDkVHR6t79+4aPHiwfHx8ytz2559/1oIFC3T48GF5enqqRYsWGj58uFq2bFnu+6Wmpmr27NnauXOnbDab4uPjNWDAAPXt27eyhwIAAOoYx6I5Mqa/J2PxXJlHPSdTWMkzNkZOlj7+fJFWNOjsHLvaOKo/9Lq6pksFUE0qfYZn+fLleuWVV7RkyRIZhqH27dsrKSlJJ0+e1IwZM/TYY4/p7Nmzpbb76KOP9NZbb+nQoUNq27atmjVrps2bN+vpp5/W2rVry3yvtWvX6umnn9amTZvUqFEjdejQQcePH9fkyZP10UcfVfZQAABAHWPq2FOKiJbSjssxcWyJMz2GtVBzP5ylr88LO+3P7tU/BrSRiWftAHVGpc/weHp6qn///ho0aJBiYmKc4xkZGZowYYL27dunDz/8UA8++KBz2datWzVv3jwFBgZq/Pjxzu127typZ555RpMnT1arVq0UEBDg3CY7O1uTJ0+Ww+HQqFGj1L17d0nSmTNn9NRTT2nevHnq3Lmz2rRpU9lDAgAAdYQpLELmUc/JMXGsM/SYRz0nhYZp9Udf6P3gbs51E3OO65GbO8srItKFFQOoapU+w9OnTx/dfffdJcKOJIWGhuquu+6SVHRmxmazOZd9++23kqThw4eX2K5Fixa67rrrlJubqyVLlpTY3+LFi5Wbm6suXbo4w44khYSE6I477pAkzZ07t7KHAwAA6pji0HP+mZ4d06ZpkldHGaaij0JhBWf1xA1J8o+KusDeANQ21dq0oFGjRpIkq9WqrKwsSVJhYaG2bt0qSerRo0epbYrH1q9fX2K8+HVZ23Tq1EkWi0VbtmxRYWFh1R0AAACoE84PPUezbXrB3kqFHhZJkp8tT0/2jVNEfLSLqwRQHao18Jw4cUKS5OHh4bw87ejRo7JarQoKClJYWFipbRo3bixJOnDgQInxgwcPSpKaNGlSahtPT08lJCTIarXq6NGjVXoMAACgbjCFReh4l+s0rv3dyvQq+lzi6bDp0S4hatwk3sXVAagulb6HpyLz58+XJHXo0EEWS9FfUdLT0yWpzLAjST4+PvL391dOTo7y8vLk6+ur3Nxc5eTkSJIaNGhQ5nYNGjTQnj17lJ6ersTExAvWNnLkyFJjXl5emjBhgiQpPDz8gvuoTp6eRf80ERH0/0dJzA1UhPmB8jA3pNTpn+vJjIY65RPiHPvnzlm6esQT8qjHXxfmBipSF+ZHtZ3h2bBhg5YsWSIPDw/ddtttzvH8/HxJReGiPN7e3iXWLf7/85ddaBsAAIBiOz6bppH7g0qEnbv2fKvex9fr9BP/lD3tuOuKA1CtquUMz+HDh/Xmm2/KMAyNGDGixBkXwzAkqcJ2j8XrVKdJkyZVuDw9Pb1G6ihPcYpOSyv/QWmon5gbqAjzA+Wpr3PDMAwd/Hq2njqToDPegc7xvzeWBl5/pxwTd8l+4ojSHr+3zOf01Af1dW7g4rh6fphMplLN0S5VlZ/hOXXqlJ5//nnl5ORo8ODBuuGGG0os9/X1lVT0wNLyFDceKH5g6fkPLi1vu+Lx8h5yCgAA6hfDMLRv1pd64kwjnfEqCjsmw6H7mpp0wxVJZXZvO/85PQDqhioNPJmZmRo/frzS09PVt29fjRgxotQ6xffGnDp1qsx95OfnKycnR/7+/s5w5OfnJz8/P0nS6dOny9yueNzV994AAADXMwxDe2ZM11PZTZ0NCkyGQ/9s4an+PVo61yP0AHVflQWevLw8vfDCCzpy5Ii6deume+65p8zL1mJjY2WxWJSZmVlm6Nm3b58kKSEhocR4cYvrvXv3ltrGZrPp4MGDslgsio2NrYrDAQAAtZRhGNr1+ed6Kq+Fsiz+kiSz4dBDSRZd261FqfV/G3qMjatqumQA1ahKAo/VatVLL72kPXv2qH379nrooYdkNpe9ay8vL7Vp00aStHr16lLLi8c6d+5cYrxTp07lbrNhwwZZrVa1adOmwmYIAACgbjMcDqV++qmetiYrx1J0dYjZcGhkK4v6dmle7nbFocd0290yXzukpsoFUAMqHXgcDodef/11bdu2TcnJyRo9erSzfV15Bg0aJEmaPXu2jh075hzfuXOnFi1aJF9fX11zzTUltunXr598fX21bt06rVmzxjl+9uxZTZs2TZI0ePDgyh4OAACopQyHXds/maZn7G2U61l0WbyH4dCY7hHq1an8sFPMFBZB2AHqoEp3aVuwYIHWrl0rSQoMDNR7771X5nojRoxQUFCQJKldu3a64YYbNH/+fD388MNq27at7Ha7Nm/eLIfDofvvv9/5oNJiAQEBuvfee/Xqq69q0qRJatWqlQIDA7Vlyxbl5ORo4MCBatu2bWUPBwAA1EKGw64tH07Tc54dlO9R9KgKT8Ouh3tGqXtT7u8F6rNKB57s7GznfxcHn7LccsstzsAjSXfeeacSExO1YMECbdmyRR4eHmrTpo1uvvlmJSUllbmPHj16aNy4cZo9e7Z27dolm82muLg4DRgwQFdffXVlDwUAANRCht2ulA8+0nNeXVToUXRpu8Ww65Ero9W1cdkPOgdQf5gMVz5sxo0dO3aM5/DALTE3UBHmB8pTV+eGYbNpw9SPNMGnmwo9LJIkL8Omx3rFqVOjUBdXVzvU1bmBquHq+VEVz+GplgePAgAAVDfDatUvUz/Wi37dZTMXfaTxNmwa2yde7RuGuLY4AG6DwAMAAGodw1qo1e99rFcCejjDjo9h1ZPXJKhNbLCLqwPgTgg8AACgVjEKC7TivY81KfAK2c0ekiRfw6qn+yUqOSbQxdUBcDcEHgAAUGsYBfla+t4nei34CjlMRWHHz7DqmesS1TKKsAOgNAIPAACoFYz8XC2ZMk1vhl4lh6noUYIBRqGeGdBUzSP8XVwdAHdF4AEAAG7PyMvVoinT9FaDq2T8L+wEGoUad30zNQ33c3F1ANwZgQcAALg1Izdb37/3uf57XtgJMgr07xuaq3EDwg6AihF4AACA2zJysjR/yhd6N6KXcyzEKNCzg1ooIdTXhZUBqC0IPAAAwC0ZWWf17ftf6v3zwk4DFejZG5MUH+ztwsoA1CYEHgAA4HaMzAx9NXWWPoq4yjkW/r+wExtE2AFw8Qg8AADArRhnTmvmB19pWuS5sBOhAo0fkqToQMIOgEtD4AEAAG7DOJ2uLz76Vl9EXukcizLl69khrRQV4OXCygDUVgQeAADgFhzpJ/TZJ9/py8iezrEYc77GD22tcD+LCysDUJsReAAAgMs5Th7TJ5/+oNmRPZxjceZ8PTu0tcIIOwAqgcADAABcynH8iKZ+vljfRnZ3jjX0KND4oW0U4stHFQCVw08RAADgMo6jhzRl+s+aH9nVOZboWaB/D22tYB8+pgCoPH6SAAAAl7AfPqB3Zi7X95FdnGNNLAUaN7SNgrw9XFgZgLqEwAMAAGqc/eBeTZ69RosiOjvHmnkVaNyQNgog7ACoQgQeAABQo2z7d+s/36zXkoiOzrEW3oV6Zkgb+XsRdgBULQIPAACoMba9qXrt2xQtC2/vHEv2KdRTQ1rLz0LYAVD1CDwAAKBGWHfu0KQF27QyvJ1zrI1voZ64sY18LWYXVgagLiPwAACAalf46zZN/CFVq8PaOMfa+Vv1xOA28vYk7ACoPgQeAABQrQq3b9aLP+7VurBWzrGOATY9Nqg1YQdAtSPwAACAalOwZaMm/HRQGxokOce6BNr1yKBW8vIg7ACofgQeAABQLfJT1un5ZceU0qClc6xbsF0PD2wli4fJhZUBqE8IPAAAoMrlbVit8atOaWtoc+fYFaEOjbq+lTzNhB0ANYfAAwAAqlTu2pV6dt1ZbQ9p6hzrFWboX/2T5UHYAVDDCDwAAKDKZK9ern9vyFZqcGPnWN9w6YHrkgg7AFyCwAMAAKpE1oqfNG5zoXYFJzrH+kWadF+/FoQdAC5D4AEAAJV2dumPGrfd0J6gBOdY/2iz7r2mucwmwg4A1yHwAACASjmz5Ac9vdNT+wPjnGM3xHror32bEXYAuByBBwAAXLbTC7/TM3t9dCAgxjl2Y0NP3dWrqUyEHQBugMADAAAuS/qCuXr6YKAOB0Q5x4YlWHTnVU0IOwDcBoEHAABcsrR53+ipow101D/COfa7xt66o2ciYQeAWyHwAACAS3Jizmw9dTJSx/3CnWO3N/XR7d0bEXYAuB0CDwAAuCiGYejE1zP15Kk4nfRt4Bz/Y3Nf3dqtkQsrA4DyEXgAAMAFGYaho7Om66mzjZTuG+oc/78kfw3v3NCFlQFAxQg8AACgQoZh6MiMz/RkTlOd9glxjv+lVYCGdox3XWEAcBEIPAAAoFyGYejQ59P0VH5LZXgHOcf/1jZIg9rFurAyALg4BB4AAFAmw+HQgc8+0VPWVjrrHegcv7d9sK5vE1PBlgDgPgg8AACgFMNh196PP9IzRjtlegVIkkyGofs6NdB1raIusDUAuA8CDwAAKMFw2LX7ww/1jLmDsr38JUlmw6H7u4TrmqRIF1cHAJeGwAMAAJwMu107P3hf4zy7KMfiJ6ko7DzUPUJ9mkdcYGsAcD8EHgAAIEkybDbtmPq+nvXuplxPX0mSh+HQyJ5RuqppmIurA4DLQ+ABAAAyrFZte/99PevbU/me3pIkT8Ou0VfEqGeTBhfYGgDcF4EHAIB6zrAWavOU9/V8wBXK9zgXdh7pFatujUIvsDUAuLcqCTx79+7V5s2btXv3bu3atUsZGRmyWCz69NNPy1z/1ltvveA+W7duraefftr5etu2bRo3bly56zdv3lzPPffcpRcPAEA9ZhQWaOO77+uFoKtU6OElSbIYdj3WO16dE4JdXB0AVF6VBJ6ZM2dq3bp1F71+nz59yl22YcMGZWVlKTk5uczlUVFRSkpKKnMcAABcPKMgX+vffV8vBvdWoYdFkuRl2DS2b0N1iCfsAKgbqiTwtGjRQomJiWratKmaNm2qv/3tbxWuf99995U5npOTo5UrV0qSevXqVeY6SUlJ5W4PAAAkx6I5MnXsKVNY+V3VjPxcrX17il4K7yebuejjgLdh0xPXJKhdbFBNlQoA1a5KAs+wYcOqYjdatWqVrFarmjdvrpgYnuAMAMClciyaI2P6ezIWz5V51HNlhh4jL1erJk/RxMhzYcfHsOnJfolqExNQ0yUDQLUyu7qA8y1btkyS1Lt3bxdXAgBA7WTq2FOKiJbSjssxcayMU2kllhu52Vo+eYpeibzWGXb8DKue6U/YAVA3uU3gSU9P16+//ioPDw9dccUV5a53/PhxffbZZ3rnnXf02WefacOGDXI4HDVYKQAA7ssUFiHzqOfKDD2OrLNaOvl9TYq6VnazhyTJ3yjUMwOaKDmSsAOgbnKbttTLli2TYRjq2LGjAgMDy10vNTVVqampJcYSEhI0atSoS7oMbuTIkaXGvLy8NGHCBElSeHj4Re+rOnh6Fv3TRETwVGuUxNxARZgfkCRFRMj+/Ns6/cQ/ZT9xRKZXn5Lj8Rc1678f6/Xo6+QwFf29M1BWvf77rkqK5p6d+oyfG6hIXZgfbhV4pPIvZ/Pz89OQIUPUvXt3Z7DZv3+/Pv/8c+3atUvjx4/Xyy+/LD8/vxqrGQAAd+UREa0G4//jDD2fv/CaJrf8nYz/hZ0gFeqNP3RTi8jy/8gIAHWByTAMo6p3euutt1b4HJ7f2rt3rx599FH5+/vr3XfflcViuej3cjgcGjdunHbs2KHbb79dw4cPv9yySzh27Jiq4Utz0YpTdFpa2gXWRH3D3EBFmB/4Lfv6lZo9f7WmNbnBORZsFOrfg1ooMdTHhZXBXfBzAxVx9fwwmUyVbmbmFvfwFJ/d6dGjxyWFHUkym80aOnSoJCklJaXKawMAoLbK/ekHvbL0cImwE1KQqWeviiDsAKg3XB54HA7HBZ+9cyHR0dGSpDNnzlRVWQAA1FqGzapjn0zVo6meWhnZzjkelXdK4zf9V/HvPl2qexsA1FUuDzxbtmxRRkaGIiIilJycfFn7yMnJkST5+PDXKgBA/WaczdCm/7yl0fb2OhAQ6xzvFOTQ+yO6KjbAs9yW1QBQF7k88BRfztarVy+ZTKbL2seaNWskSY0bN66yugAAqG0ce1P1zTuf6d/h/ZVt8XeOD42RXruztxrEx5XbshoA6iqXBp6CggKtXbtW0oUvZ1u4cKGysrJKjBmGoYULF2revHkymUzq379/tdUKAIA7y1u2SK99s0EfxPVztp32clj1UGsf/eWaJHmai/6oWNFzegCgLqqSttQbNmzQrFmzSozZbDaNHTvW+frmm29Wp06dSqzzyy+/KD8/X02bNlVcXFyF7/H1119r6tSpio+Pd3aLOHjwoE6ePCmTyaQ777xTTZo0qYrDAQCg1jBsNp2cPk0TsuK1N7Kjczy84Kwe7dtQzZvEltqmOPQ4Jo6V0o7L2LhKpmuH1GTZAFBjqiTwZGZmateuXSXGDMMoMZaZmVlqu/MvZ7uQwYMHKyUlRYcPH9aWLVtkt9sVGhqqXr16aeDAgWrWrFkljwIAgNrFyDyjLR98qJeDrlJmYIBzvLXOaMz1jRUaG1XutsWhx9i4SmbCDoA6rFqew1MX8BweuCvmBirC/Kg/HPt3ad6MBfog5hrZzR7O8UEJPvrLlY2cl7AVY26gPMwNVMTV86MqnsNTJWd4AABAzclfuUTvrj6qH+Ouc45ZDLvu6Rala1uEu7AyAHA/BB4AAGoJw25X+pef6cWMSO2K6uwcDzNb9ci1TdUyws+F1QGAeyLwAABQCxhZmdr+wVS9FHCFzgQFOceT/B16dECyQn35lQ4AZeGnIwAAbs44tE/ffzZHU2Kvlc187lf3gHgv/fWqxrJ4XN5z7ACgPiDwAADgxgrWLNX7y/fp+/gBzjFPw6G/do7Q9ckRLqwMAGoHAg8AAG7IcNh1etYXeimtgX6N7u4cDzHZ9Mh1jdUq0t+F1QFA7UHgAQDAzRg5WUqdOlUv+nfX6eAQ53hzP7seG9BSYX4W1xUHALUMgQcAADdiHN6vRZ99rXdi+slqPhds+sVadE/vFvLyMLuwOgCofQg8AAC4Ceu6lfrgp1TNi7veOeZhOHRXx3Dd0CpCJhPNCQDgUhF4AABwMcNh15mvv9QrRwO1NaanczzIZNPD/RqrbTT36wDA5SLwAADgQkZutnZ/MFUTfLoqPTTUOd7Ex67Hr2+pCH/u1wGAyiDwAADgIsaxQ/rp45maHNNPhR5ezvE+0Z66r08LeXtyvw4AVBaBBwAAF7BtXK2PF23TN/EDnWNmw9D/tQvV0LZR3K8DAFWEwAMAQA0yHA5lfvulJh70VUrslc7xANk05ppG6hAb6MLqAKDuIfAAAFBDjLxc7f3wfb1o6aQTDcKc44nedj1+fQtFBXhVsDUA4HIQeAAAqAHG8SNa/vF0vRl9rQrOu1/nykgPPXB1C/lwvw4AVAsCDwAA1cyW8os++36jZsXd4BwzGYZGtAnW8PYx3K8DANWIwAMAQDUxDEPZc2dp0j4PbYjr7Rz3k02j+jZSl3ju1wGA6kbgAQCgGhj5eTr40ft6wdxex8IinOPxXnaNHdBCsUHcrwMANYHAAwBAFTNOHtPqjz7X65HXKM/TxznePdysh65pLj+LhwurA4D6hcADAEAVsm/ZoOnz12p6/A0lxn+fFKhbO8XKzP06AFCjCDwAAFQBwzCU893Xen2XQ2vj+zrHfWXXv3rFq3tCsOuKA4B6jMADAEAlGQX5OvzxVE0wWutweJRzPNZi09gBzRUf7O3C6gCgfiPwAABQCUbacf3y4TS9GtlPuZ6+zvEuDcwa2S9Z/l7crwMArkTgAQDgMjm2p2jmtyv0WdwNMkznHhx6S/MA/aFrHPfrAIAbIPAAAHCJDMNQ7g/f6s0d+VoVf41z3Ed2PXBlnK5MDHFdcQCAEgg8AABcAqOwQMc+eV8TbMk6ENHCOR7ladPj/ZspMdSngq0BADWNwAMAwEUyTqVpwwcfaVL4Ncr28XOOtw8xacy1yQr05n4dAHA3BB4AAC6C49et+vrrn/RJ/A1ynHe/zrCm/vpTt3h5mLlfBwDcEYEHAIAKGIah/MXz9dbmTC1reK1z3Et23d8zVr2bhLqwOgDAhRB4AAAoh2Et1IlpH2hCQTPti2zqHI/wKLpfp0kD7tcBAHdH4AEAoAzG6XRt/uBDvdKgrzIDA5zjbYOkh69LUpAPv0IBoDbgpzUAAL/h2Lldc2cv0gfxA+UwnWtEMLixn/7co6E8uV8HAGoNAg8AAOfJX7JA/92QriUN+zvHLIZd93aPVr/mYS6sDABwOQg8AABIMqxWpX3+oV7MaaTd0V2c42EeNj12XVM1D/N1YXUAgMtF4AEA1HvGmdPaNnWqXg7tozNBgc7x5EBDj16XpBBffl0CQG3FT3AAQL1m7PlV381YoPfjr5fNfO7X4sBGPrqrZyNZPLhfBwBqMwIPAKDeKlj6g6asPa6FCdc7xzwNh/7eNUr9W3K/DgDUBQQeAEC9Y9isOjX9E710NlapMd2c46Fmmx7p10TJkX4urA4AUJUIPACAWs+xaI5MHXvKFBZR4XrGqTQZqxYrddchvRTcS6eDg53LWvgberR/S4X5Waq7XABADTK7ugAAACrDsWiOjOnvyTFxrIxTaeWuZ5xKk2PCw1q0JlVPhF+v097nwk6/eG89fyNhBwDqIgIPAKBWM3XsKUVES2nHyw09xqk0FT43WlPCr9JbSbc6mxN4GA79rWOY7u+dKIsHvxIBoC7ipzsAoFYzhUXIPOq5ckOP4+QxnX7uUT3T9DZ9F3+lczzYbNO/r2ukQa0iZDLRiQ0A6ioCDwCg1isv9DgO7tHul1/Qw63v1vaQJs71m/o5NHFIS7WJ8ndh1QCAmkDTAgBAnVAcehwTxxaFnsf/pp/D2+nttn9Voce5e3P6xnrpH70S5e3J3/wAoD7gpz0AoM4whUXI9ODTKvDw1vtNBun1Vr93hh2z4dBd7Rvoob6NCTsAUI9UyRmevXv3avPmzdq9e7d27dqljIwMWSwWffrpp2WuP2PGDM2cObPc/Q0dOlR//OMfy1yWmpqq2bNna+fOnbLZbIqPj9eAAQPUt2/fqjgUAEAtZqRu0YaZ3+jdLg/phO+5B4cGyqrR/RqrQ0yAC6sDALhClQSemTNnat26dZe8XcuWLRUdHV1qvEmTJmWsLa1du1aTJk2SYRhKTk5WYGCgtm7dqsmTJ+vAgQP6v//7v0uuAQBQ+xk5Wcr48lN9cDpISxN/V2JZYvZRPXJ0vmKuf0QSgQcA6psqCTwtWrRQYmKimjZtqqZNm+pvf/vbRW3Xr1+/iz4zk52drcmTJ8vhcGjUqFHq3r27JOnMmTN66qmnNG/ePHXu3Flt2rS53MMAANQyhmHIsWapfvxxrT6Ku0bZ5zUhMBsODW3krVsXzJT3ycNyTBwr86jnLvhwUgBA3VIlFzEPGzZMt956qzp37qyQkJCq2GUpixcvVm5urrp06eIMO5IUEhKiO+64Q5I0d+7canlvAID7MdKO6/B/JumpdVl6K/FGZVvOhZ3mvnZNGtREd/ZqKt+R4y74nB4AQN1Va+7aXL9+vSSpR48epZZ16tRJFotFW7ZsUWFhYU2XBgCoQYbdroIFX2n6ezP1r5D+2hrazLnMx16gvyb56sVhrdQ41EfShZ/TAwCo21zalnrr1q3av3+/CgsLFRYWpo4dO5Z7/87BgwcllX1/j6enpxISErRnzx4dPXpUiYmJ1Vk2AMBFjP27tH36l3o79AodSkgusazr2V3625CuiowvfW/ob1tWGxtXyXTtkJoqGwDgQi4NPEuXLi3xevr06erevbvuu+8++fj4OMdzc3OVk5MjSWrQoEGZ+2rQoIH27Nmj9PT0iwo8I0eOLDXm5eWlCRMmSJLCw8Mv9jCqhadn0T9NRATXmqMk5gYqUlfnhyMvVyc+fU/v7szX9/HDSywLsxi6r8FpXfuHgfKMjCl/JxERsj//tvJX/yz/G2+r5ordT12dG6g85gYqUhfmh0sCT3R0tEaMGKGOHTsqPDxcOTk52rFjh6ZNm6Y1a9bI4XBozJgxzvXz8/Od/+3t7V3mPovHz18XAFD75f2yXN9/PkdTovooIzbIOW6SoWHJ4bq3bwsFeF/crzOPiOh6GXYAoD5zSeDp3bt3idc+Pj666qqr1Lp1a40ePVq//PKLUlNT1bJly2qrYdKkSRUuT09Pl2EY1fb+F1KcotPSuM4cJTE3UJG6ND+MsxlKm/6J3s2L0y8JN5ZYluDj0D96Jyo5wk95mRnKc1GNtUldmhuoWswNVMTV88NkMikmpoKz9xfBrZoWhIaGOttUp6SkOMfPv7ytoKCgzG2Lx89fFwBQ+xgOh6w/L9C3//lQ93v30S/hrZ3LLHLoj21CNWlYspIj/FxYJQCgtnDpPTxlKU5wGRkZzjE/Pz/5+fkpNzdXp0+flp9f6V9yp0+fluT6e28AAJfPOHpQe774TP/17aRdiQNLLGsbYtY/ejVRbJCXi6oDANRGbhd4srOzJZU+U9OoUSPt2LFDe/fuVXx8fIllNptNBw8elMViUWxsbI3VCgCoGoa1UPnzZmn6r5n6Jv5GOUwezmUBZrv+3DVW/ZqGyGQyubBKAEBt5FaXtBmGoV9++UVS6fbTnTp1kiStXr261HYbNmyQ1WpVmzZt5OXFX/4AoDYxUrdqwysT9WB6gr5q2KdE2OkT563JN7XUtc1CCTsAgMtS44EnMzNTP//8s6xWa4nx/Px8TZkyRbt27VJISIi6detWYnm/fv3k6+urdevWac2aNc7xs2fPatq0aZKkwYMHV/8BAACqhJGTrYyP39Grc1P074SbdMI3zLksyuLQM9c01Mi+jRXs43YXIwAAapEq+S2yYcMGzZo1q8SYzWbT2LFjna9vvvlmderUSfn5+Xrrrbc0depUxcfHKywsTLm5udq3b5+ysrLk7++vkSNHlmo/HRAQoHvvvVevvvqqJk2apFatWikwMFBbtmxRTk6OBg4cqLZt21bF4QAAqpFhGHKsXaYlC1frg/h+yo72dy4zy9DQFkH6fccYeXu61UUIAIBaqkoCT2Zmpnbt2lVizDCMEmOZmZmSpMDAQA0dOlS7du3S8ePHtX//fpnNZkVGRqpPnz4aPHhwuQ8X7dGjh8aNG6fZs2dr165dstlsiouL04ABA3T11VdXxaEAAKqRkX5Chz//RO+YWmpL4yElljUPkP7Rq7GaNKDbJgCg6pgMVz5sxo0dO3aM5/DALTE3UBF3nR+G3a7CRXP09frD+rJhX1nNFucyH5Ndd3SI0g1JYfIwc59OdXHXuQHXY26gIq6eH1XxHB4ujAYAVCvjwG7tmP6lJof01KFGySWWdY3w1N+vbKoIf0s5WwMAUDkEHgBAtTDy85T99XRNO2jo+9ihMkzn7skJ9bDrrz0b6oqEQLqvAQCqFYEHAFDljM2/aOW3P+q9mKt1Oi64xLLrE/00omucArw8ytkaAICqQ+ABAFQZ42yG0qZ/oil5sVqbOKzEsoY+Dt3XK1HJkX6uKQ4AUC8ReAAAlWY4HLItW6gFy7fp0/hrlOd/rtOaRQ7d2iZMN7WJlMWDy9cAADWLwAMAqBTj2CHt/fwzve3bQbsSbyixrE2IWf/o1URxQV4uqg4AUN8ReAAAl8WwWpU/f6Zm7Dirb+IGy24+d09OgMmuP3eLVb+mITQlAAC4FIEHAHDJjJ1btXHWt3onopeON+xYYlnvWG/d1bOhQnz4FQMAcD1+GwEALpqRk60zsz7VB2kB+jnhphLLIi0O3XtVgjrFBrioOgAASiPwAAAuyDAMOX5ZriULV+nDuH7KivZ3LjPL0JAWQfp9xxj5eJor2AsAADWPwAMAqJBx6qSOfD5N76i5NicOKbGsmb90X+/GatLAp5ytAQBwLQIPAKBMht0u66Jv9fX6Q/oy/joVelicy3xk1x87RmlQUpg8zDQlAAC4LwIPAKAU48Ae7Zg+Q28H99DBRkkllnUJ99Q9VzVVhL+lnK0BAHAfBB4AgJORn6ecb6Zr2gGHFsQOlWE6d09OqIddf+3ZUFckBNJqGgBQaxB4AACSJGPLOq36dpGmRPfV6biQEssGJPrpT13jFODlUfbGAAC4KQIPANRzRmaG0qdP05ScaK1pNKzEsngfh+7rlahWkX6uKQ4AgEoi8ABAPWU4HLItX6Tvl23VtPhrlOd3rtOapxy6pXUD3dw2UhYPWk0DAGovAg8A1EPGscPa98Wnetu7g3Ym3lBiWetgs/7Ru7Hig7xdVB0AAFWHwAMA9YhhtSp//izN2HFG38QNlt187p6cAJNdd3aN1bXNQmhKAACoMwg8AFBLORbNkaljT5nCIipczziVJmPjKpkSmmrTrG/034jeOt6wQ4l1esV46+6eDRXiy68FAEDdwm82AKiFHIvmyJj+nozFc2Ue9Vy5occ4lSbHy48pMzNbHzYdrJ8ShpdYHmmx696rGqlTbEBNlA0AQI3jTlQAqIVMHXtKEdFS2nE5Jo6VcSqt1DqO9JOyPzdKP3nG6/5uY/RTdBfnMrMMDWseqDeHJxN2AAB1GoEHAGohU1iEzKOeKzf0OPbs0LHnn9C4xr/TG8m3K8vi71zW1F+aOLCx/twtTj6e/BoAANRtXNIGALVUcehxTBzrDD22Z/+jrB/nadaafZrR4X4Velic6/vIrj90iNLg5DB5mGlKAACoH/jTHgDUYs4zPeFFZ3pWjvmX/rEvRNOa3FAi7HQJ99SbQ1toaOtwwg4AoF7hDA8A1GKGYUhHD+hASEN9FdFXyyI7yDCd+1tWiIddd/eI11WNgmg1DQColwg8AFALGQ67tGGVUn/8SbO8k7Q2/pZS6/SPNuv/rmquAG+PMvYAAED9QOABgFrEsNnkWP2ztixbrVkBbbU5/nel1onPPaF7U2cp2Stf5g7PSd4VP6cHAIC6jMADALWAUVggx/JFWrdyk2Y16KzUhJtLrdMo76TuvLad+iY0U8YTnzobGVT0nB4AAOo6Ag8AuDEjP1e2Jd9p5fpdmhXRTQcalw46zTMP6HcZG9Xtb39RVHJzSSrVvY3QAwCorwg8AOCGjOxMFS6ap5+2HNJX0VfoWOPWpdZpl7VfN+/5Xm0sOfL4TaApq2U1oQcAUB8ReADAjRhnTiv/hzlauPOUvo65UqeadCy1TrdIi2627VHznyZLEdHlBpnfhh5j4yqZrh1SE4cBAIDbIPAAgBsw0o4r+/s5mn8gT3Njr1Rm44ASy80y1CvOVzd3iFGjEG9JTeWw5MvUsWeFZ22KQ4+xcZXMhB0AQD1E4AEAFzKOHVLG/Dmae9Kk72J7KjfRt8RyTzl0TWKghreLUkygV4llFxtgTGERnNkBANRbBB4AcAHjwG6d/G6uvjkboIUxvVSYUDLMeMuhAc1DNKxNhML8LC6qEgCA2o/AAwA1yNi5TYe/n6+vCqL0c9S1sgWW/DHsb7JrUHK4bkwOU5APP6IBAKgsfpsCQDUzDEPatkF7flik2aZGWhUxSA6TucQ6wR52DW0TqYEtG8jP4uGiSgEAqHsIPABQTQyHQ9q4Wjt+XKpZ3i20LmpoqXUiLA7d1C5a1zYLkbenuYy9AACAyiDwAEAVM2w2OdYs1aalazQ7sJ22xg0rtU6ct0PDO8aqb+NgeZpNNV8kAAD1BIEHAKqIYS2UffkirV21WbNCO2l3wk2l1mniZ+h3nePUIz5QHgQdAACqHYEHACrJyM+V7afvtWz9Ls2O6KZDjYaVWic5yKRbOsWpU6y/TCaCDgAANYXAAwCXycjJUsGieVq85bC+ju6pE4nJpdbp2MBDt3SOU+tIPxdUCAAACDwAcImMsxnK/WGOvt95RnNieiqjcYcSy00y1CPKS7/rGKdmYT6uKRIAAEgi8ADARTPSTyjz+281/0C+5sb2VHaif4nlZhnqE++rmzvEqGGwt4uqBAAA5yPwAMAFGMcO6/R3czTnpIcWxPRQfqOSYcYih65tHKib2kUpKsDLRVUCAICyEHgAoBzGwT06/t08fX02QD/G9JK1oaXEch+TXdc3D9WwNpEK9eXHKQAA7qhKfkPv3btXmzdv1u7du7Vr1y5lZGTIYrHo008/LbWuw+FQamqq1q9fr+3bt+vkyZPKzc1VWFiY2rZtq2HDhikyMrLUdtu2bdO4cePKraF58+Z67rnnquJwANRzxq7tOvj9As0uiNTSqH5yBHiUWB5gtuvG5HANSg5XoLdHOXsBAADuoEoCz8yZM7Vu3bqLWvfkyZN6+umnJUkNGjRQixYtZDabtXv3bi1atEgrVqzQY489pqSkpDK3j4qKKnNZVFTU5R8AgHrPMAxp20btWvSjZqmR1oRfL8NkLrFOqIddQ9tG6foWDeRrMZezJwAA4E6qJPC0aNFCiYmJatq0qZo2baq//e1vFa7fvn173XTTTWrVqpVzzGq1asqUKfrpp5/0xhtv6I033pCnZ+nykpKSdN9991VF2QAgw+GQNq3W1h+Xa6ZPS22KuLHUOlEWu4Z3iNE1TUPk5UHQAQCgNqmSwDNs2LCLXjc6Olpjx44tNW6xWHT33Xdr7dq1Sk9P186dO0sEIgCoSobdLseapdqw7BfNCmirHbFDSq3T0MehmzvGqndisDzMPCwUAIDayK3usvXy8lJMTIz27Nmj06dPu7ocAHWQYS2UbcVirV61WbNCO2tffOmg08zf0C2d49UtPkBmE0EHAIDazK0Cj8PhUHp6uiQpJCSkzHWOHz+uzz77TFlZWQoMDFRSUpI6dOggs5nLTIC6zrFojkwde8oUFlHhesapNBkbV8l87bkwY+Tnyfrz9/p5/R59FdFVRxJKB502wSbd0jle7aP9ZCLoAABQJ7hV4FmxYoXOnj2roKAgtWzZssx1UlNTlZqaWmIsISFBo0aNUkxMzEW/18iRI0uNeXl5acKECZKk8PDwS6i86hXfvxQRUfEHO9Q/9XVu5Hw7XVnT35P5p+/UYPx/5BERXeZ69rTjOv3qU3KcOCK/gAD5Xj1QGXNnac7a3foqsrvSG5VuetIz2kd39m6ptrFB1X0Y1a6+zg9cGHMD5WFuoCJ1YX64TeBJT0/Xhx9+KEm69dZbZbGUfN6Fn5+fhgwZou7duzuDzf79+/X5559r165dGj9+vF5++WX5+fnVdOkAaoBPjz7Knful7CeO6PQT/ywz9NjTjuv0E/+U/cQRmSOidfbwEX381CuaE91DZxOalVjXJENXJwTqT1c1V4vIgJo8FAAAUINMhmEYVb3T4sBS1nN4ypKfn69nnnlGe/fuVdeuXTVmzJiLfi+Hw6Fx48Zpx44duv322zV8+PDLLbuEY8eOqRq+NBetOEWnpaW5rAa4p/o8N4xTaXJMHCulHZciomUe9Zzz8rbzl2UGhGteZBfNj+mhHEvJP4J4yKG+Df10c4dYxQV5ueIwqlV9nh+oGHMD5WFuoCKunh8mk+mSruIqi8tvfLHZbJo4caL27t2rpKQkPfjgg5e0vdls1tChQyVJKSkp1VEiADdhCouQedRzUkS0lHZcjolji+7XOZUmx4uP6NTZXH3QdLD+3vEhfZlwTYmw4yWHBjfx1zvDmuuB3ol1MuwAAIDSXHpJm8Ph0JtvvqmUlBQ1atRIjzzyiLy8Lv1DSHR00WUtZ86cqeIKAbib4tBTfDbH8fjfdNw7RF817KMlbbvIZi75Y83PZNfAFqEa0iZSIT5ucxUvAACoIS797f/ee+9p1apViomJ0RNPPCF/f//L2k9OTo4kycfHpyrLA+DG7J2uVMovW7UkurNWRbSTw1TyhHWQ2a4bW4XrhuRwBXh5uKhKAADgai4LPJ999pkWLVqk8PBwPfnkkwoODr7sfa1Zs0aS1Lhx46oqD4CbMXKy5Vi3XKkbt2tpQbBWRLZXZrvupdYL87BrWLso9W/RQD6eLr9qFwAAuJhLAs/cuXP19ddfKyQkRE8++eRFtYBeuHChevToocDAQOeYYRhatGiR5s2bJ5PJpP79+1dn2QBqmGEtlDb/ooNrN2jpGQ8ti2ivExE3lLluTG66bspYr6v//Ht5Rbi2rTwAAHAfVRJ4NmzYoFmzZpUYs9lsGjt2rPP1zTffrE6dOmn//v365JNPJEmRkZGaPXt2mfvs16+fkpLOPS/j66+/1tSpUxUfH+/sFnHw4EGdPHlSJpNJd955p5o0aVIVhwPAhQyHXUrdqvQ1q7T0aKGWNWitfUHXSWU8IsfTYVPnBmb1TQxSl09flUfaMenVbTLO694GAADqtyoJPJmZmdq1a1eJMcMwSoxlZmZKKrrfprjd886dO7Vz584y99m6desSgWfw4MFKSUnR4cOHtWXLFtntdoWGhqpXr14aOHCgmjVrVuZ+ALg/wzCkQ3uVvXq5Vuw5pWWBLbQtpLeMRqUvSTMZhlqf2aPeuXt1xZ9uU2B0VNE+Ro0/18hg4tgSLasBAED9VS3P4akLeA4P3FVdmhtG2nEVrFmmX7Yf0jKvhlofllyqy1qxJl6F6nVgpa7at1xhwX5lBpqKntNTX9Sl+YGqxdxAeZgbqIir50dVPIeHHq0AapSRlSnbL8u1JWWnltnDtDqijXLj2pa5bpTFrt5Nw9SnWajifvlexg/zKwwyv21ZbWxcJdO1Q6r7kAAAgBsj8ACodkZBgRyb1mjP+hQtzfLV8oh2yogs+567ILNdVzUKUt8W4WoR5iOTyVS04NohckgydexZ4Vmb4tBjbFwlM2EHAIB6j8ADoFoYdrv062YdXbNWS0/atSysjY4EXyeV0YHeR3Z1j/FR36QotYv2l6fZVOY+LzbAmMIiOLMDAAAkEXgAVCHDMKT9u5WxeoWWH8jU0uAk7QrqIzUsva6HHOoY6qE+ydHq1jCQZ+YAAIBqQeABUGnGyaPKWbVMq3ce1zKfRG0O7SZHw7IDTLK/Q32So3Rlo2AF+fAjCAAAVC8+bQC4LEbmGRWuXa4NW/ZqqaK0LqyVCmPbl7lugpdNfVpGqFeTUEUFeNVwpQAAoD4j8AC4aEZ+nuwbV2v7hu1amhegVeFtlR1Z9jOwwj1s6t0kRH2ahysx1KeGKwUAAChC4AFQIcNmk3Zs0r41G/TzaZOWh7VRekg/KaT0ugEmm66M91efpCglR/jKbCq7+QAAAEBNIfAAKMUwDGlvqk6sXqWlh/O0LCRZBwN6S3Gl1/WSXd0iLOrTKkYdYwJk8SDkAAAA90HgAeBkHDuss6uXa8XudC31b6pfg3tK8aXXMxuG2gcb6t0qVj0SAuVn8aj5YgEAAC4CgQeo54wzp5W3drnWbjukpR6x2hTaVvbYsgNMCx+beidHq1fjEIX48uMDAAC4Pz6xAPWQkZcr24ZV2rRpp5YWhGhteCvlR7Yoc904T6t6Nw9Tn+bhigmkwxoAAKhdCDyAG3MsmiNTx54yhUVUuJ5xKk3GxlUyXzuk/HVsVhlbNujXdZu19IynVoS3VWbI1WWuG2q2qldCoPokRatpA2+ZaD4AAABqKQIP4KYci+bImP6ejMVzZR71XLmhxziVJsfEsVLacTmkEqHHcDikPb/q4Oq1WnrcqmWhrXTCr5fkV3o/frKpZ7S3+rSOVZtIP3mYCTkAAKD2I/AAbsrUsaeMxXOLgszEsWWGnvPDjiKiZerYs2j8yEGlr16hpfvOallAc+0L7CHFln4PT8OhLg1M6t06Vl3jA+XlYa6JQwMAAKgxBB7ATZnCImQe9dy5szf/Cz2KKAo9pcLOXaOUvXaFVvx6XMss8doW0l5GTOkAYzIMtQmwq3erGF2RGKIALzqsAQCAuovAA7ixskKP/fm3Jck5VhgYqnWR7bXs2y1aH5YkW2RSmftq4lWo3i0j1btZmML8LDV5GAAAAC5D4AHc3G9DT9o9t8huSFuDE7Ws5S1aHdFGuZ6+ZW4bZS5UnyYh6p0UpYbB3jVcOQAAgOsReAA3ZzgcUm62srr00dZftiulQXP9EtZKGd5BZa4fZLKqV5yf+rSOVYswHzqsAQCAeo3AA7gh48wpFW7bpNTUg9p02qEU/4baG9hRjjady1zfRzb1CPdU7zbx6hDjT4c1AACA/yHwAG7AKCiQsXOrDmzbqU3HcrTZI1zbQpqowDtOiil7Gw+HXR0z96p3r3bq0aqhvD3psAYAAPBbBB7ABQyHQzq0T6e2btGm/ae0udBPKcFNdca7oxRZ/nYhhdnqHO2nbs0ilfTxvxV4fL90JLqoe9sFHk4KAABQHxF4gBpinE5T7tYUbdt1RJvOSpsDGumQf5IUUv42XoZdrTL3q33aDrXXaTW+7yFFJbeWJJ1s8GSpltXlPZwUAACgviLwANXEyM+TPXWrdm3fo00n8rTZEqWdQQmy+cRLPmVvYzIMNfUqUPu4ILUP91bLj8fLcvKIFBFdKtCU95weQg8AAMA5BB6gihgOu4z9e3Rs63ZtOnRGKbZAbQlpolzPDlIFGSTKVKD2EV5q3yxG7WIDFeRd9CBQx6I5MsoJO8V+G3qMjatkunZINR0hAABA7UPgASrBOHVSZzenaMue49qU7aGUwESl+VR8mZq/rGobaKhD0yh1SAhVTKBXmeuZrx0ihyRTx54VnrUpDj3GxlUyE3YAAABKIPAAl8DIy1Xhji3asWOfNqVbtdk7RnsD4mT4NpTKfvanPA2HWnoXqkNCiNo3jVKzBj4X3Tb6YgOMKSyCMzsAAABlIPAAFTDsdjn27dL+rb9q05EsbXYEa3twogo9Okjh5W/X0JynDlG+at88Vm2iA+RroWU0AACAKxB4gN8w0o4rbfMWbdp3UptzvLQ5KFFnvSq+TC1UBWofbFL75jHqkBCqBr58awEAALgDPpWh3jNys5W7bYu2pB5Symm7UnxidcS/oVTBZWrehk1tfK1qnximDk2jlBDsJZPp4i5TAwAAQM0h8KDeMWw22fbu1K6tu7TpeI5S1EC7AhvKbm4nhZW9jdlwqKlnvjrE+Kt9y1glRQTI4kHAAQAAcHcEHrgdx6I5F+xMJknGqbSL6kxmGIaME0d1ZMt2bdp/Sin5PtoalKg8z5ZScPnbRStPHRp4qH2LOLWLD1HA/9pFAwAAoPYg8MCtOBbNkTH9PRmL51b4EE3jVNq5B26qdDczIydLZ7duUcrOI0o5I6X4xSvdp6Hk11DyK/u9A4xCtfW3qUPjSHVsFqmogLLbRQMAAKD2IPDArZg69pSxeG5RkJk4tszQc37YUUR00TY2qwp2p2r71r1KOZmvFHOY9gXGSeYGUoOy38vTsCvJkqcO8UHq0CJeTcJ8L7pdNAAAAGoHAg/cSvFDNJ1nb34Tekqc2QkN176ug5Qye4k2F/prR1CCrOYWFV6mlqgctYvwUoeW8WodFywfT9pFAwAA1GUEHrid8kKPkZOlk29NVIpXQ6W06qctIU2VmRtQdIlaOZepNTDy1T7QoQ7NotS+SZRCaRcNAABQr/DpD27JFBYh3TVKJ6f8R7+aQrXj/W+1ObS5jrW7v8LtfBxWtfHOV4eEELVPaqiGwd60iwYAAKjHCDxwC4bdLuuhfdr9636lHsv8//buPTiq+u7j+Htz3VwJgVzBkBBgFxrkphgfHastFRHrIFSknVrp1Dpe2ipgO8yDgzLCY0d9oJcZnM5jUQpSawUdBUw7FssYK3eQQJLNJpCbhEtCwsKGXDa7zx90U0KyK5Czm+z6ec3sOJ7Lb3+/mS8n57PnnN+h/GIUNnMmzROf8rtfhMfNmAgnkzPMTLLmYMlM1nTRIiIiItJNgUcGhKfVyVmbjfJjDZQ3tWNzJ1KVkI0rYiTEcOnjQ1bXBSYNNTF5XDYT89JJjNF00SIiIiLSNwUeCTiPx0PX6ZMcL6vEVn+W8vMmbNHDOR2XCuRDgv/9s1obsZyvYfy5aiY12UhPNvudslpERERExEuBRwzncXXiqKzCZq+j/Ewrto5Y7PHZtEemgSkNkn3vG+N2MSbCiSXRg+XIDsZ9WUJKShIRS1YB+Jy9TURERESkLwo80m9djnN8WVpBWc0Zys91YTMN4cv4dOAGiOXSx4fhXa1YzB1YM5Owjr2BvPQkoloae7xn5/Jg42/KahERERGRKynwyDXxuN1c/LKeivLjlDc4KL8YTUVsOs7oJCDJ7+1pkZ4uRnvOYxkSgTUnDUt+NulJPdPQlS8VvTLQfNV7ekRERERELqfAI36529o4VVFJWVUD5WcvTS5QG5eO25QBkRmQ6Hvf5K6LWKMvYk2LxzJmBGNHDiP2K1706Tn4uc+w43Vl6PEc/BzTjPv7O1QRERERCUMKPNJDe1MjVUcrKfuy+d+TCwzjXEwScAOYfe9n8njIcTuwJrixjhiK1ZJDVkrcNb8DJ2LG/bgB05Rb/V618YYez8HPiVDYEREREREfFHi+xjzuLpqqqimz12M77aS808xxczquiBQgxe/tafFd7VgiLmBJjWH86EzGjs4iIcaYcrraAGMalqYrOyIiIiLilyFnqMeOHePw4cNUVlZit9tpbm4mOjqat956y+9+O3fupKioiPr6eqKiohg3bhxz587FYrH43Mdms7FlyxYqKipwuVyMHDmSmTNncueddxoxlAHh/viDr7yiAZeeb+nPFY1Op5PqUvtlkwuk0BibAmRANJc+PoxwncNidmHNTMRqGcUNGUOIuMarNyIiIiIiwWZI4Hn33XfZt2/fNe2zfv16tm3bRkxMDDfeeCOdnZ0cPnyYL774gsWLFzN9+vRe++zZs4fVq1fj8XgYP348SUlJHDlyhLVr11JTU8MjjzxixHCCyv3xB3j+8jqeHVv9Pnx/+cP8br76KojH48FxooHy8hrKTzoovxhFZUwaHZFmvur2tFh3B2M9DizJEVhHpWGxjGJIvJ83gYqIiIiIDFKGBJ5x48aRm5tLfn4++fn5PPbYY363P3LkCNu2bSMpKYmVK1eSlZUFQEVFBS+88AJr165lwoQJJCb+54n4CxcusHbtWtxuN0uWLOGWW24BoKWlheXLl7Nt2zamTZtGQUGBEUMKGtOUW/Hs2Op3xrErZy4zTbm1VztdnR3U2Y5Rfuwk5U0dlLsTaTCnAkMvfeJ89yGt8zzWmItYhsdjHTuSvFEZREXo6o2IiIiIhD5DAs+cOXOuafsPP/wQgLlz53aHHbgUnL7zne/w0Ucf8cknn/Dd7363e92OHTtobW3lpptu6g47ACkpKfzwhz/k1VdfZevWraEXeHxMs0zapdDja5pmZ3MLFaXHKK8/S/kFExVRw2iNMgOZ4OdiTJS7i9FdzVgTPVhHDMUyPpfhKX6mWhMRERERCWFBn7Sgo6ODI0eOAFBYWNhrfWFhIR999BH79+/vEXj279/vc5+pU6cSHR1NSUkJHR0dxMSE1u1XfYWerv95DQD3/y7Dc+YkDRljsd18L7Zth7B1xlEbm4rHFA/E+709LcXlxBJxAevQGMbnZzJ6zA3ERmuuChERERH5egj6me+JEyfo7OwkOTmZYcOG9Vqfl5cHQE1NTY/ltbW1AIwePbrXPlFRUeTk5FBVVcWJEyfIzc01vuMBdnnoaW9q4p+/ehZbcg62jJnYxo7CEZMIF4DIf3/6EOFxM6qzGYu5k/FZSVituWRkWK55amgRERERkXAR9MDT2NgI0GfYATCbzSQkJOB0Orl48SJxcXG0trbidDoBSE1N7XO/1NRUqqqqaGxsvKrAs3jx4l7LYmJi+PWvfw3A8OHDr2Y4xkpLo+t/XuOFV9/iH1m9J224UoLrIuNN5ykYFsONY0dQMHk8iQl+HtaRsBAVdemfbVqa/1n95OtJ9SG+qDbEF9WG+BMO9RH0wNPW1gbg97az2NhYnE4nbW1txMXFde/jXedrn8vbD2WW87V9Bp4bOs4y3tzJxOwhTJqYT97YXCIjIgaghyIiIiIioSHogcfj8QD4vc3Ku00grV692u/6xsbGoPTjct4JCqxON+audsY66rA4arC4mrA8/DBD8v6rx/Znm5qC2j8ZHLy/sJw5c2aAeyKDkepDfFFtiC+qDfFnoOvDZDL1mOTsegQ98MTFXbrlqr293ec2HR0dwKXb2y7/r3e/+Pj4Xvt427t821By+WxsI4Znsv37VqKjCzjz309cmqHt/+rx+HlPj4iIiIiI9Bb0+6G8z8Y0+bg60dbWhtPpJCEhoTscxcfHd4ecs2fP9rmfd/mAPHvTT1dOPR357CrM2SOJ/Pc01KRlds/e5mnSry8iIiIiIlcr6IEnOzub6OhoHA5Hn6Hn+PHjAOTk5PRYPmrUKACOHTvWax+Xy0VtbS3R0dFkZ2cHoNeB4+s9O17e2dsUekRERERErl3QA09MTEz3y0F37drVa7132bRp03osnzp1qs99Dhw4QGdnJwUFBSH3Dh7Pwc99hh2vK0OP5+DnA9BTEREREZHQMyBTfM2ePRuALVu20NDQ0L28oqKCjz/+mLi4OL71rW/12Ofb3/42cXFx7Nu3j927d3cvP3fuHBs3bgTgvvvuC0LvjRUx435MDz3qM+x4eUOP6aFHiZhxfxB7KCIiIiISukweA6YiO3DgAJs3b+7+f7vdjslkYsyYMd3L5s2b132VBuDNN99k+/btxMbGMnHiRLq6ujh8+DBut5tFixZRWFjY63t27drFmjVrAJgwYQJJSUmUlJTgdDqZNWsWP/7xj/s7lG4NDQ1Bn6XtcgM9I4YMXqoN8Uf1Ib6oNsQX1Yb4M9D1MWhmaXM4HNjt9h7LPB5Pj2UOh6PH+oULF5Kbm0tRURElJSVERkZSUFDAvHnzsFqtfX5PYWEhK1asYMuWLdjtdlwuFyNGjGDmzJncddddRgxFRERERETCiCFXeMKRrvDIYKXaEH9UH+KLakN8UW2IPwNdH0Zc4RmQZ3hERERERESCQYFHRERERETClgKPiIiIiIiELQUeEREREREJWwo8IiIiIiISthR4REREREQkbCnwiIiIiIhI2FLgERERERGRsKXAIyIiIiIiYUuBR0REREREwpYCj4iIiIiIhC0FHhERERERCVsKPCIiIiIiEraiBroDg5XJZBroLgCDpx8y+Kg2xB/Vh/ii2hBfVBviz0DVhxHfa/J4PB4D+iIiIiIiIjLo6JY2EREREREJWwo8g9TSpUtZunTpQHdDBiHVhvij+hBfVBvii2pD/AmH+tAzPINUR0fHQHdBBinVhvij+hBfVBvii2pD/AmH+tAVHhERERERCVsKPCIiIiIiErYUeEREREREJGwp8IiIiIiISNjSe3hERERERCRs6QqPiIiIiIiELQUeEREREREJWwo8IiIiIiISthR4REREREQkbCnwiIiIiIhI2FLgERERERGRsKXAIyIiIiIiYUuBR0REREREwlbUQHfg66Kjo4P333+fzz77jMbGRhITE5k0aRIPPfQQw4YNu6a2nE4nf/3rX9mzZw8tLS2kpKRw8803M3/+fBISEgI0AgkUI2rD6XRy8OBB9u/fT3V1NY2NjZhMJkaOHMntt9/O3XffTVSU/rmHIiOPHZdraGjg2WefpbOzk0mTJrFs2TIDey3BYHRtnDx5kvfff5+SkhJaWlowm81kZWUxffp07r///gCMQALFyNo4dOgQ27dvp6qqitbWVhISEhgzZgyzZ89m4sSJARqBBMKxY8c4fPgwlZWV2O12mpubiY6O5q233rqu9kLpfNTk8Xg8A92JcNfR0cGLL76IzWZj6NChWK1Wzpw5Q2VlJcnJyaxcuZLMzMyrauv8+fM899xzNDQ0kJGRwejRo6mvr6euro7MzExWrVpFUlJSgEckRjGqNt5++222bNmCyWQiLy+PzMxMHA4HNpuNzs5OrFYry5YtIzY2NgijEqMYeey40ooVKygtLcXj8SjwhCCja2PPnj389re/xeVykZubS1ZWFhcuXKC2tpbY2Fh+//vfB3A0YiQja2Pr1q386U9/wmQyYbFYSE1N5dSpU1RVVQHw6KOPcvfddwdyOGKgl19+mX379vVYdr2BJ9TOR/WTbxC899572Gw2xo0bx3PPPYfZbAb+cyB57bXXWLFixVW1tX79ehoaGpg+fTqLFi0iMjISgHXr1lFUVMT69ev52c9+FrCxiLGMqg2z2cwDDzzAzJkzSU1N7V7e0NDAiy++SHl5OZs3b+YHP/hBwMYixjPy2HG5HTt2cPToUWbMmMHHH39sdLclCIysjerqan7zm98QFxfHL3/5S6xWa/c6t9vN8ePHAzIGCQyjasPhcLBp0yaioqJYvnx5j7rYtWsXa9asYcOGDdxxxx3d3yGD27hx48jNzSU/P5/8/Hwee+yx624r1M5H9QxPgLlcLoqKigD4yU9+0uOgcN999zFq1CjKyso4duzYV7bV0tLCp59+SmRkJI8++mh3cQE8/PDDJCcnU1xcTEtLi+HjEOMZWRtz5szh+9//fo+wA5CVldUdcj777DMDey+BZmR9XO7cuXNs2LCBiRMncttttxnaZwkOo2vjjTfewOVy8eSTT/Y4qQWIiIggPz/fuM5LQBlZG3a7HZfLRUFBQa+6KCwsJCcnh/b2durr640dhATMnDlzmD9/PtOmTSMlJeW62wnF81EFngArLy/H6XSSkZFBXl5er/W33HILQK9LjH05ePAgHo+HCRMm9CrU6Ohopk2bhtvt5tChQ0Z0XQLMyNrwJzc3F4Dm5uZ+tSPBFaj6eOONN+jo6OCnP/2pIf2U4DOyNurr6ykrKyMrK4tp06YZ3lcJLiNrIzo6+qq+MzEx8do6KSEvFM9HFXgCrKamBqDPAw/A6NGje2zXn7a8y6urq6+1mzIAjKwNf06dOgXQr19zJPgCUR8HDhzgX//6Fw888MB1P/sjA8/I2jhy5AgAN954Ix0dHfzzn/9k3bp1rFu3jn/84x+0trYa1GsJBiNrIz8/n/j4eI4cOUJ5eXmPdbt376a2thaLxaJjyddQKJ6P6hmeAGtsbATwOSuK9xYk73ZX09aVty15eb/jatqSgWdkbfizfft2AG666aZ+tSPBZXR9tLW18cc//pHs7GzmzJljSB9lYBhZG3V1dQDExMTwq1/9ihMnTvRYv2nTJpYsWcKECRP602UJEiNrIyEhgccff5zf/e53PP/8892TFpw+fZqqqiomT57Mk08+aVznJWSE4vmoAk+AtbW1AficHct7f613u/605V3e3t5+zf2U4DOyNnz5+9//TklJCQkJCTrJDTFG18fbb7/NmTNnWL58uaYoD3FG1obT6QQu/TCSkJDAs88+S0FBAS0tLbz77rsUFxfzyiuvsHr1aoYOHWrQCCRQjD5uFBYWkpiYyJo1a3pc5RkyZAjf+MY3BtUsXBI8oXg+qr96AfZVs35fy6zg3m1NJlO/+iSDg5G10ZfS0lLefPNNTCYTTzzxhM9fYmRwMrI+qqqqKCoq4o477qCgoKC/XZMBZmRtuN1uALq6uvj5z3/OpEmTAIiPj+cXv/gFDQ0NVFVV8be//Y0FCxZcf6clKIz+u/Lhhx+ycePG7nerpKenc/r0af7yl7+wceNG7HY7S5Ys6U+XJQSF4vmonuEJsLi4OMB3yvUuv5opHb1t+fplxtuW3rUSGoysjSvV1NTwyiuv4HK5WLhwIdOnT7/+jsqAMKo+urq6+MMf/kB8fDw/+tGPjO2kDAgjjx3ebVJTU7vDzuXuuusuAI4ePXpdfZXgMrI2SktL2bBhA7m5uSxevJicnBzMZjM5OTksWbKEvLw8du/ezRdffGHcACQkhOL5qK7wBNjw4cMBaGpq6nP92bNne2x3NW1597mS9zuupi0ZeEbWxuVOnjzJqlWrcDqdPPjgg8yaNat/HZUBYVR9NDU1UV1dTUpKCqtXr+6xzns7U2VlJS+88AJms5mlS5f2t+sSYEYeO9LT0wFIS0vrc713ucPhuOZ+SvAZWRs7d+4ELs3sFhHR8/fxiIgIpk+fzvHjxzl69GifYVnCVyiejyrwBNioUaMAfL64zTsXvne7/rTlXX41bcnAM7I2vM6ePcvKlStpaWnh3nvv5cEHH+x/R2VAGF0fLS0tPt+J4HQ6KS0tJT4+/to7KkFnZG14p62/cOFCn+vPnz8PXN+VZgk+I2vDezLr/TX/St7lvmpHwlcono8q8ASY1WolPj6eU6dOcfz48V5T+O3evRuAqVOnfmVbkydPxmQyUVZWxrlz5xgyZEj3us7OTvbv34/JZGLKlCnGDkICwsjagEt/dFatWsXp06e58847eeSRRwzvswSPUfWRnp7OO++80+e6o0ePsmLFCiZNmsSyZcuM6bgEnJHHjokTJxIbG8vJkydpbGzs9YtsaWkp4Hv6WRlcjKwN7zlGVVVVn+u9y71XCeXrIxTPR/UMT4BFRUVxzz33ALBu3boe9ztu3bqVmpoarFYrY8aM6V5eVFTEM888w6ZNm3q0NXToUG677TZcLhevv/46XV1d3es2btyIw+Hg9ttv1/tWQoSRtdHe3s5LL71EXV0dt956K48//nhIPUwovRlZHxJejKyN2NhYZs2aRVdXF6+//nqPtg4dOsTOnTsxmUzMmDEjwKMSIxhZG95nP4uLi3u9qHTv3r0UFxdjMpn0jGgYC6fzUV3hCYK5c+dSUlKCzWbj6aefxmq10tjYiN1uJykpqdc89g6HgxMnTtDc3NyrrYULF2K329m9ezfPPPMM+fn51NXVUVdXR0ZGhn7VDzFG1caf//xn7HY7ERERREZG8tprr/X5fU899VTAxiLGM/LYIeHFyNr43ve+R1lZGQcOHODpp59mzJgxOBwOKioq8Hg8LFiwoMcJsgxuRtXGzTffTGFhIbt27eLll18mPz+ftLQ0zpw50311Z8GCBWRnZwdtbNI/Bw4cYPPmzT2WuVyuHlf4582b130FMJzORxV4giAmJobnn3+e9957j+LiYvbu3UtCQgLf/OY3eeihh67poa7k5GReeukl3nnnHfbu3cuePXsYMmQI99xzD/PnzycxMTGAIxGjGVUb3ofP3W43xcXFPrdT4AktRh47JLwYWRvetj744AM+/fRTDh06RHR0NAUFBcyePfuqb6uVwcGo2jCZTCxatIhPPvmEnTt3UltbS3V1NfHx8UyZMoVZs2YxefLkwA5GDOVwOLDb7T2WeTyeHsuudoKSUDsfNXn6+7IPERERERGRQUrP8IiIiIiISNhS4BERERERkbClwCMiIiIiImFLgUdERERERMKWAo+IiIiIiIQtBR4REREREQlbCjwiIiIiIhK2FHhERERERCRsKfCIiIiIiEjYUuAREREREZGwpcAjIiIiIiJhS4FHRERERETClgKPiIiIiIiELQUeEREREREJWwo8IiIiIiISthR4REREREQkbCnwiIiIiIhI2Pp/axh1emxcqIQAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Conditions\n", "alpha = 0.1\n", "ti = 100\n", "ts = 300\n", "\n", "t_target = 1.0\n", "dt = 0.01\n", "\n", "# Make grid\n", "nx = 11\n", "x = np.linspace(0, 1, nx)\n", "dx = np.diff(x)[0]\n", "\n", "# Const\n", "adx2 = alpha / dx**2\n", "\n", "# Solution array\n", "u = np.ones_like(x)*ti\n", "du = np.zeros_like(u)\n", "\n", "# Calculation\n", "t = 0\n", "while abs(t - t_target) > 1e-10:\n", " # Adjust time step to reach target time\n", " dt = min(dt, t_target - t)\n", "\n", " # Apply bc\n", " u[0] = ti\n", " u[-1] = ts\n", " \n", " # Spatial discretization\n", " cs(nx, u, du, adx2*dt)\n", " \n", " # Update solution and time\n", " u += du\n", " t += dt\n", "\n", "# Exact solution\n", "dts = ts - ti\n", "u_exact = ti + dts*x + np.sum([\n", " 2*dts*(-1)**n / (n*np.pi) *np.exp(-n**2*np.pi**2*alpha*t)*np.sin(n*np.pi*x) \n", " for n in range(1, 10)\n", "], axis=0)\n", "\n", "# Visualization\n", "plt.plot(x, u, marker='x')\n", "plt.plot(x, u_exact)\n", "plt.legend(['Computation', 'Exact'])\n", "\n", "# Compute Error\n", "err = np.linalg.norm(u[1:-1] - u_exact[1:-1])\n", "err = u - u_exact\n", "err = np.sqrt(np.sum(err**2) / nx)\n", "print('Error: {:.5f}'.format(err))" ] }, { "cell_type": "markdown", "id": "a625f806", "metadata": {}, "source": [ "### Accuracy, Stability Convergence\n", "\n", "#### 정확도\n", "FTCS 의 Truncation Error 는 $O((\\Delta t), (\\Delta x)^2)$ 이다.\n", "즉 시간 차분에 대해서는 1차 정확도, 공간 차분에 대해서 2차 정확도를 갖는다. \n", "\n", "Exact solution 과의 오차는 다음과 같이 계산한다.\n", "\n", "$$\n", "L_2 Error = \\sqrt{\\frac{1}{N} \\sum (T-T_{exact})^2}\n", "$$\n", "\n", "선형 편미분 방정식이므로, Linear Wave 방정식과 같이 Consistency, stability 분석이 가능하다.\n", "\n", "또한 이를 통해 수치 해석 결과가 이론 해로 수렴 (Convergence) 하는지 확인할 수 있다.\n", "\n", "#### von Neumann Stability\n", "\n", "FDE (Finite Difference Equation) 의 완전해를 $D$ 라 했을 때 수치해 $N$은 다음과 같이 나타낼 수 있다.\n", "\n", "$$\n", "N = D + \\epsilon.\n", "$$\n", "\n", "여기서 $\\epsilon$ 은 round-off 에 의한 error 이다.\n", "\n", "이를 차분식에 적용하면, $D$ 은 차분식을 만족하므로 사라지고, $\\epsilon$ 만 남는다.\n", "\n", "이 오차를 다음과 같은 형태로 표현하자.\n", "\n", "$$\n", "\\epsilon_j^n = \\sigma^n e^{ikx_j}\n", "$$\n", "\n", "이를 차분식에 적용하면\n", "\n", "$$\n", "\\frac{\\sigma^{n+1} e^{ikx_j} - \\sigma^n e^{ikx_j}}{\\Delta t} \n", "=\n", "\\alpha \\frac {\\sigma^n e^{ikx_{j+1}} -2 \\sigma^n e^{ikx_j} + \\sigma^n e^{ikx_{j-1}}}{\\Delta x^2}\n", "$$\n", "\n", "정리하면\n", "\n", "$$\n", "\\sigma = 1 + \\frac{\\alpha \\Delta t}{\\Delta x^2} (2 \\cos (k\\Delta x) - 2)\n", "$$\n", "\n", "$|\\sigma| < 1$ 을 만족하기 위해서는\n", "\n", "$$\n", "\\frac{\\alpha \\Delta t}{\\Delta x^2} (2 \\cos (k\\Delta x) - 2) \\geq -2\n", "$$\n", "\n", "즉 \n", "\n", "$$\n", "\\Delta t \\leq \\frac{\\Delta x^2}{\\alpha (1 - \\cos (k\\Delta x))} \\leq \\frac{\\Delta x^2}{2 \\alpha}\n", "$$\n", "\n", "### Convergence\n", "\n", "- Consistency와 Stability를 만족하므로 Lax Equivalance Theorem에 의해 수치해가 이론해에 수렴함을 알 수 있다." ] }, { "cell_type": "markdown", "id": "f62a089d", "metadata": {}, "source": [ "### Du Fort-Frankel Method\n", "FTCS에서 시간에 대해서도 Central Difference를 하면\n", "\n", "$$\n", "\\frac{T_i^{n+1}- T_i^{n-1}}{2 \\Delta t} \n", "=\n", "\\alpha \\frac {T_{i+1}^n -2 T_i^n + T_{i-1}^n}{\\Delta x^2}\n", "$$\n", "\n", "이 기법은 Unconditionally unstable 하다.\n", "\n", "수치 안정성을 확보하기 위해 우변 항을 변화하면 Du Fort-Frankel 기법이다.\n", "\n", "$$\n", "\\frac{T_i^{n+1}- T_i^{n-1}}{2 \\Delta t} \n", "=\n", "\\alpha \\frac {T_{i+1}^n -2 \\frac{T_i^{n+1} + T_i^{n-1}}{2} + T_{i-1}^n}{\\Delta x^2}\n", "$$\n", "\n", "이 기법은 Unconditionally stable 하다.\n", "\n", "Talyor expansion을 이용해 Truncation Error을 확인하면\n", "\n", "$$\n", "\\frac{\\partial T}{\\partial t} - \\alpha \\frac{\\partial^2 T}{\\partial x^2}\n", "=\n", "-\\frac{\\Delta t^2}{6}\\frac{\\partial^3 T}{\\partial t^3}\n", "+\\frac{\\alpha \\Delta x^2}{12}\\frac{\\partial^4 T}{\\partial t^4}\n", "+\\frac{\\alpha \\Delta t^2}{\\Delta x^2}\\frac{\\partial^4 2}{\\partial t^2}\n", "+...\n", "$$\n", "\n", "$\\Delta t \\rightarrow 0$, $\\Delta x \\rightarrow 0$ 으로 점근하더라도 Truncation Error는 사라지지 않는다.\n", "\n", "즉 Consistency를 만족하지 않는다." ] }, { "cell_type": "markdown", "id": "c274181e", "metadata": {}, "source": [ "## 실습\n", "1. FTCS 기법을 완성하시오.\n", "\n", "2. FTCS 기법에 대해서 $\\Delta t$ 를 바꿔가면서 von Neumann Stability Analysis 결과를 확인하시오.\n", "\n", "3. $\\Delta t=10^{-5}$ 로 매우 작게 정하자. 이때 $\\Delta x=0.1, 0.05, 0.025, 0.0125$ 로 바꿔가면서 오차를 구하시오. 공간에 대해 2차 정확도를 만족하는지 확인하시오.\n", "\n", "4. $\\Delta x=0.05$ ($n_x=21$) 로 정하자. 이때 $\\Delta t=0.01, 0.005, 0.0025, 0.00125$ 로 바꿔가면서 오차를 구하시오. 시간에 대해 1차 정확도를 만족하는지 확인하시오." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.2" } }, "nbformat": 4, "nbformat_minor": 5 }